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I. Summary of Research for the Third Quarter 

Research activities during this quarter for comets have focused upon 
completing two major papers describing the gas coma of comets. Research 
activities for the hydrogen distribution in the Saturn system have concen- 
trated upon assessing exciting new data uncovered in Voyager 1 and Voyager 2 
data tapes. Originally scheduled work to perform preliminary Titan torus 
model calculations was therefore postponed. 

1 . Cometary Coma 

Two major papers were completed this quarter and submitted to The Astro- 
physical Journal for publication. Both papers are attached in the appendix of 
this report. 

The first paper is entitled "Monte Carlo Particle Trajectory Models for 
Neutral Cometary Gases. I. Models and Equations." This paper documents a new 
class of Monte Carlo particle trajectory models based upon physical processes 
of atoms and molecules in the coma. The model provides a new level of realism 
in predicting the spatial distribution of observed chemical species in come- 
tary comae and in exploring the basic physics of the transition zone between 
true fluid- flow and free molecular flow. 

The second paper is entitled "Monte Carlo Particle Trajectory Models for 
Neutral Gases. II. Spatial Morphology of the Lyman-a Coma." This paper 
successfully applies the general model of the former paper to simulate the 
hydrogen Lyman-a image for Comet Kohoutek and includes near perihelion one 
image representing an extreme case for collisional thermalization. The 
success of the model in this application is noteworthy. 

2. New Observational Data for Saturn-System Hydrogen 

A significant effort of the past two years has been extended in acquiring 
additional data to characterize the hydrogen distribution in the Saturn 
system. The main avenue to pursue this goal has been the collaboration effort 
established with D.E. Shemansky, who is a member of the ultraviolet spectro- 
meter instrument team (UVS) of the Voyager spacecraft. In this period, a 
search of the magnetic tapes containing all pertinent UVS data ofj Voyager 1 
and Voyager 2 was undertaken. The search has uncovered a significant amount 
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of new data that had not previously been reduced. In fact, this new data set 
is approximately 2.5 times larger than the old data set that was already 
reduced and published by the UVS team. A summary of the older published data 
and newly reduced data is given in Table 1. 

There are a number of interesting aspects of the new data in Table 1 that 
deserve discussion. The most spectacular is that from the Voyager 1 post- 
encounter data an image of the H Lyman-a emission has been constructed at a 
viewing angle of —25° to the orbit plane of Titan. The image indicates that 
the brightness distribution is not cylindrically symmetric as expected for a 
Titan torus alone. The brightness distribution between Saturn and Titan's 
orbit has a minimum at pre-dawn and a maximum in the vicinity of the dusk 
terminator line suggesting an asymmetric dayside hydrogen source for a large 
electroglow- driven planetary corona. In the less corona- dominated portions of 
the circumplanetary hydrogen, the presence of the Titan hydrogen torus can be 
distinctly identified. In addition to these data, the newly reduced Voyager 2 
pre- encounter data provide the best one -dimensional data in the orbit plane. 
Combining the earlier published data and new data in the analysis of the H 
distribution with an appropriately developed model for the planetary corona 
and the Titan hydrogen torus model should provide a most interesting study. 
This study will be initiated here and continued in our renewal proposal to be 
submimtted to NASA early in the next quarter. 

II. Program of Research for the Next Quarter 

Research activities in the next quarter will focus upon (1) further 
assessing the new data for the hydrogen distribution in the Saturn system, and 
(2) integrating new plasma data for the planetary magnetosphere currently 
undergoing reduction by E.C. Sittler into the Titan torus model to improve the 
lifetime description of hydrogen. 


Table 1 

Voyager Spacecraft Data Available for Research on the Atomic 
Hydrogen Distribution in the Saturn System 


Voyager 1 pre -encounter data 

Published Data: 14 days of one -dimensional scan data nearly 

parallel to the satellite orbit plane have been reported in the VI 
30 day science report (Science 212, 201-211, 1981). These data were 
obtained at a range of - 2 x 10 7 km from Jupiter. 

New Data: Additional data of similar nature is available but the 

quality is uncertain. The reduction process should be repeated to 
include the additional data. 

Voyager 1 post-encounter data 

Published Data: none 

New Data: Observations were obtained in mosaic scans across the 

system. All the data have been reduced for the first time in the 
past few months of 1987. The Voyager 1 post encounter trajectory is 
out of the solar system plane, and the data are obtained at an angle 
of ~25° to the orbit plane of Titan. An image of the system in H 
Lyman-a emission has been constructed but is not yet published. The 
data in this image were obtained over the period 1980 DOY 324-343 
(6-25 days post encounter) at spacecraft planet ranges of 0.83 x 10 7 
- 0.33 x 10 8 km. The integration time of the data set is 126 
hours. Titan moves - 1.25 orbits during the 20 day observing 
period. 

Vovager 2 pre- encounter data 

Published Data: Observations of one -dimensional scan data 

approximately normal to the satellite orbit plane were published in 
the V2 30 day science report (Science 215, 548-553, 1982). 

New Data: All usable observations for one -dimensional scan data 

approximately parallel to the satellite orbit plan have been reduced 
for the first time. The data were obtained in 1981 DOY 180-186, at 
ranges of 0.54 x 10 8 - 0.49 x 10 8 km. The integration time for this 
data set is 88 hours. These observations provide the best one- 
dimensional data in the orbit plane of Titan that were obtained by 
the Voyager spacecrafts . 

Vovager 2 post -encounter data 

Post -encounter sequences were not obtained because of the lockup of 
the scan platform drive. 
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Abstract 


The Monte Carlo particle trajectory models for the distributions of 
neutral gases in cometary atmospheres, developed in original form several 
years ago have been greatly extended and generalized. In a Monte Carlo 
particle -trajectory model (MCPTM) the spatial distribution of a neutral 
cometary species is determined by calculating the explicit trajectories of 
many individual particles following the chain of parent vaporization, coma 
outflow, photodissociations, collisions, and decay through the appropriate 
number of generations. The original inner coma model has been generalized and 
now includes: (1) solar radiation pressure, (2) isotropic ejection of 
daughters owing to excess photolysis energy, (3) radially variable parent 
outflow speed and gas temperature, (4) time (heliocentric distance) dependent 
source rate and lifetimes, (5) heliocentric velocity dependent lifetimes and 
excitation, (6) multiple collisions with the outflowing gas which includes 
both the effects on the observed species as well as the energy transfer 
(photochemical heating) delivered to the coma. An extended coma MCPTM has 
also been developed for application to the spatially extended atomic clouds 
(e.g. H, C, and 0), which uses the generalized inner coma MCPTM for the source 
region and calculates explicitly the atom trajectories in three dimensions, 
including the relative heliocentric orbits of the nucleus and cloud atoms and 
heliocentric velocity dependent radiation pressure where important. In. this 
paper the mathematical derivations of the various methods employed in the 
MCPTM are presented in detail, and its application to the calculation of the 
photochemical heating of the inner coma through the partial thermalization of 
cometary hydrogen atoms produced by the photodissociation of water is 


discussed. 


I. Introduction 


Models of the neutral comae of comets have generally focused on two 
extreme limits of actual conditions which occur in cometary comae. These are: 
(1) exospheric models which describe the distributions of particles in the 
outer coma where collisions between particles are rare, and (2) fluid based 
models which describe the bulk properties of various fluids (composed of 
molecules, atoms, dust etc.) that are to one extent or another collisionally 
coupled to each other or are at least self-coupled, and which are thus used to 
investigate such complex issues as gas -phase chemistry and/or dynamics. 

The first models were simple exospheric models used to describe the 
observed distributions of gas and dust in the comae of comets. The fountain 
model of Eddington (1910) described the space and column density distribution 
for particles emitted isotropically from a point source (the cometary nucleus) 
moving under the influence of a constant radiation solar radiation pressure 
force. Further enhancements to the basic fountain model were added and are . 
documented in the work of Wallace and Miller (1957). Perhaps the most notable 
contribution was that of Haser (1957,1966) whose models for radial outflow of 
a daughter species produced by the photodissociation of an unknown parent take 
into account either production and decay scale lengths or one scale length and 
radiation pressure acceleration. The two scale length model is still largely 
used today, despite many shortcomings, to analyze the observations of neutral 
radicals and atoms. 

Models introduced by Keller and Thomas (1975) and Keller and Meier (1976) 
were isotropic point source models for application to the large spatially 
extended hydrogen coma as observed above the earth's atmosphere in the light 
of resonantly scattered solar Lyman-cr. These models used the calculation of 
the syndyname on the sky plane to approximate the apparent distortion 



introduced by the fact that the H atoms were actually following heliocentric 
orbits under the influence of solar gravity and solar radiation pressure. 

Meier et al. (1976) were able to use this model to analyze two-dimensional 
sky-plane images of comet Kohoutek (1973 XII). They required the use of the 
weighted superposition of three Maxellian distributions (with most probable 
speeds of A, 8 and 20 km sec’^) to describe the effective outflow speed 
distribution of H atoms from the inner coma region ( < 10^ km) . Since the 
speed distribution resulting from the excess energy of the photodissociation 
of water (which based on circumstantial evidence was suspected to be the 
ultimate source of observed H atoms and OH radicals) was calculated to have a 
speed distribution sharply peaked at 8 and 20 km sec'^, Meier et al. suggested 
that the thermal izat ion of H atoms in the inner coma by collisions with heavy 
parent gas was responsible for the low speed atoms. 

An important departure for the exospheric models took- place with the 
introduction of models which take into account the (non-radial) isotropic 
ejection of daughter radicals (or atoms) upon photodissociation of the 
parent. In Independent work using very different methods Combi and Delsemme 
(1980a) and Festou (1981a) developed new models which take into account this 
isotropic ejection. The main impetus behind this type of work was the fact 
that parent (source) scale lengths determined by applying Haser * s model to the 
observed radial brightness distributions of cometary radicals were found to be 
generally smaller than values calculated using estimated cometary gas outflow 
speeds and photodissociation lifetimes computed from the solar uv flux and 
measured photoabsorption cross sections. 

Combi and Delsemme (1980a) developed two approaches to this problem. One 
is the average random walk correction to the interpretation of Haser* s 
model. The other is the explicit simulation of the isotropic ejection by 


Monte Carlo methods. Their models were applied to the observed spatial 
distributions of CN in comets (Combi and Delsemme 1980b, Combi 1980). Festou 
(1981a) developed the vectorial model in which the isotropic ejection of 
photodissociated OH radicals from the outflowing ^0 molecules is explicitly 
solved for and numerically integrated. Both studies concluded that Haser's 
model though useful in describing the apparent distribution of observed 
radicals, underestimates the true scale lengths of true parent molecules. 

More recently, Kitamura et al. (1985) have added a multiple collision 
algorithm to the basic Monte Carlo model of Combi and Delsemme (1980a). They 
have used their new model to begin the study of the collisional thermalization 
of cometary hydrogen atoms in the inner coma and conclude that collisions are 
important for the combination of small heliocentric distance and large gas 
production rate. Later in this paper we will return to the discussion of 
their work since we have considerably extended and improved their initial 
approach. It is important to note here that the addition of collisions to a 
basically exospheric model represents a first step toward modeling the 
transition zone between collision-dominated fluid flow and essentially 
collisionless free molecular motion. Bockelee-Morvan and Gerard (1984) and 
Schloerb and Gerard (1985) have also adopted the Monte Carlo method rather 
than the vectorial model, since they need to describe both the space density 
and the velocity distributions of OH necessary for interpreting 18-cm radio 
observations in comets. 

The fluid-based models have been applied to the collision dominated inner 
coma. Originally, two types of modeling in this area had been developed. 

These were the pure hydrodynamic model and the Lagrangian gas -phase chemical 
model. In the pure hydrodynamic models the dynamics (flow speed and 
temperature) of the coma was the object of study. The paper by Mendis and Ip 


(1975) gives an excellent review and a complete documentation of the early 
work in this area. The Lagrangian gas -phase chemical models assumed the coma 
•could be described as a single fluid parcel traveling away from the nucleus in 
which non- equilibrium gas phase chemistry (in addition to photochemistry) 
could occur. Important examples of such models were developed by Giguere and 
Huebner (1978), Mitchell et al. (1981) and Cochran (1985). Typically these 
models contain the possible inclusion of many species (>100 including parent 
molecules and their fragment radicals, atoms and ions) and a large network of 
chemical and photochemical reactions (>1000 ). More recently, dynamics and 
chemistry have been combined in chemical gas-dynamic and dusty -gas -dynamic 
(Marconi and Mendis 1982, 1983, 1984; Huebner and Keady 1983; Ip 1983; 
Crovisier 1984; Gombosi et al. 1985). An excellent review article by Gombosi, 
Nagy and Cravens (1986) presents a detailed discussion of the methods and 
procedures involved in such models. 

What typically distinguishes the models of different investigators from 
one another is the extent to which different processes are included and by 
what methods they are included, for example: IR, visible and UV radiative 
transfer, gas -phase chemistry, radiative cooling, photochemical heating, 
multiple fluids for various classes of gas species, electrons, ions and dust, 
and finally full time -dependent dynamic equations. Even two dimensional 
axisymmetric dusty-gas models have been developed to study the evolution of 
gas and dust jets (Kitamura 1986, and Komle and Ip 1986). 

As mentioned earlier the fluid based models are applicable deep in the 
collision zone where the gas can reach thermal equilibrium and can be 
described at any point in space as one or more fluids (being coupled or not) 
each having in general a bulk velocity and one or more temperatures (i.e., 
kinetic and internal states) . Departures from local radiative equilibrium 


have been treated by iterative methods in most of the gas -dynamic models (see 
references above). However, when the collisional mean free path for molecules 
becomes too large the applicability of a hydrodynamic model breaks down even 
for the flow mechanics. The gas in these models is described in terms of a 
phase space distribution consisting of a bulk motion and a kinetic tempera- 
ture. A sufficient collision rate is required in order to constantly 
re-equilibrate a Maxwell-Boltzmann distribution and eventually transfer this 
internal energy to bulk flow through adiabatic expansion. The free -flow 
exospheric models are applicable in the outer regions where essentially no 
collisions between molecules occur. A characteristic distance separating 
these two regions is typically defined such that the distance to the nucleus 
is equal to the local mean free path for intermolecular collisions. This 
distance called the collision zone radius, is the direct analogue to the 
exobase level in planetary atmospheres (Chamberlain 1963) . 

Figure 1 shows the collision zone radii for a number of representative 
comets as a function heliocentric distance , as estimated from published gas 
production rates and/or visual light curves. It is seen that the collision 
zone radius can vary from values of less than 100 km (as will be expected for 
much of the time CRAF will be in the vicinity of comet P/Tempel 2) to values 
greater than 100,000 km for comets Kohoutek and West near their perihelion 
times. For much of the time extensive ground-based observations were made of 
comet P/Halley (0.6 - 2 AU) the collision zone radius was of order one to a 
few tens of thousands of km. Also plotted in Figure 1 are several values of 
the percentage of water photodissociated as a function of distance from the 
center of the nucleus and of the heliocentric distance. What is critical here 
is that the higher the percentage of water dissociated within the collision 
zone, the higher the heating rate and the larger the extent of hydrogen atom 


thermalization will occur. The simple relation between the collision zone 
radius and the water photodissociation scale length is complicated by the 
large mass ratio between the hot freshly produced H atoms and the outflowing 
heavier gas atoms and molecules so many collisions are required to thermalize 
the H atoms completely. Therefore, comparison of collision zone radii and 
percentage of water photodissociation only provides an upper limit to the 
heating and thermalization. 

In reality, there is no clear separation between the inner (collisional) 
and outer (collisionless) regimes, and a large transition region exists at 
distances of order of the collision zone radius. For the less productive 
periodic comets (P/Encke and P/Tempel 2) , it appears fairly clear that an 
exospheric (free -flow) model is in fact applicable to the generally observable 
comae and that the fluid dynamic models are only applicable out to a few tens 
of km at best. However for the large productive comets, fluid models should 
be applicable in a much larger inner coma (perhaps a few thousand km) , but 
typical distances of importance to observations are made right in the middle 
of the transition region. Furthermore, even in the case of the planetary 
exobase where , owing to the presence of gravity and the exponential decrease 
in density with increasing height, the transition from a fluid to a 
collisionless corona is much more spatially confined, the applicability of a 
thermalized Maxwell -Boltzmann distribution is questionable (Fahr and Shizgal 
1983) . In the case of a cometary atmosphere the density decreases 
approximately with an inverse square law so the transition is much more 
gradual . 

The Monte Carlo particle trajectory model (MCPTM) allows us to bridge the 
gap between the regimes of the purely fluid to the purely free flow models and 
to investigate the conditions in this all-important transition region. In 


fact, the MCPTM provides the only reasonable way to describe accurately the 
appropriate space -time -velocity distributions for molecules, radicals and 
atoms even in otherwise purely free flow models, owing to such complications 
as solar radiation pressure, orbital motion effects, time dependences, noniso- 
tropic production, and correct inner coma outflow conditions. 

In the remainder of this paper we will describe the mathematical 
formalism used in the application of the principle of Monte Carlo to a wide 
variety of processes in the cometary atmosphere. Section II describes the 
basic procedures used in the inner coma model. Section III describes the 
geometrical and orbital mechanics procedures for the extended coma model. 
Section IV describes a particular application which couples a steady- state 
inner coma MCPTM with a simple gas -dynamic model to study the photochemical 
heating of the coma by the thermalization of photodissociated hydrogen 
atoms. Section V presents a summarized discussion of the implications of the 
MCPTM, its capabilities, as well as future work to be pursued in this area. 

The companion paper to this one (Combi and Smyth 1987b hereafter called Paper 
2) describes in detail the coupled aspects of the thermalization of hydrogen 
atoms as it relates to the spatial morphology of the observed two-dimensional 
distribution of the Lyman -a coma and to the dynamics of the inner coma. Some 
of the results of using the MCPTM have been described previously (Combi and 
Smyth 1985, Combi, Stewart and Smyth 1986, Combi 1987, and Combi and Smyth 
1987a), however, most of the mathematical and numerical details concerning the 
MCPTM appear only in this and the companion paper in this issue. 



II. Monte Carlo Methods for the Inner Coma Source Model 


The original application of the Monte Carlo method to describe the 
spatial distribution and kinematics of neutral gases in cometary atmospheres 
was discussed by Combi and Delsemme (1980a). Further refinements and addi- 
tions have been discussed by Bockelee-Morvan and Gerard (1983) and Kitamura et 
al. (1985). In this section of the paper we present a number of fundamental 
model advances, some of which are in the interest of model efficiency or 
improved statistics (such as the use of skewed probability density distribu- 
tions and appropriate weights, or forced dissociations), while others are to 
add generalizations or new physical processes (such as time dependent life- 
times, Maxwell -Boltzmann distributions, or multiple particle collisions). 

In all applications the Monte Carlo method is used by computing the actual 
trajectories of many (-10^) radicals or atoms, making sure that all of the 
important variables (time, space, velocity etc. ) are appropriately spanned, 
and then calculating space and column densities, and radiative emission rates 
by explicitly counting the total weighted numbers of particles (radicals or 
atoms) within regions of specified volume and area, respectively. Similarly, 
doppler profiles can be computed by also counting particles in specified line- 
of -sight velocity bins. 

In all situations, the principle of Monte Carlo is used to span the pos- 
sible states a particle may occupy relative to some variable or condition. 
These variables or conditions can be related to the weight associated with a 
particle, owing to production rate variations, or source and decay lifetimes, 
or they can relate to the trajectory calculations, owing to isotropic redirec- 
tion upon dissociation or collision, or a thermal ized speed component. Many of 
the derivations in this paper are direct analogues or extensions of processes 
outlined in a monograph by Cashwell and Everett (1959), which was written to 



address a number of problems in photon and neutron scattering in bulk material 
having complex geometrical configurations. In order to use the Monte Carlo 
method, a state or property of a particle, defined symbolically as x' , must 
be chare terized in terms of a probability density distribution function, 
p(x')dx'. A probability distribution function can then be defined such that 

P(x) - J X p(x')dx', on some interval from x- to x^. If the probability 
X 1 x 2 

density is properly normalized on this interval, then x^J p(x')dx'— 1, and 

values, x^, of the distribution can be spanned according to the relation 

x. 

PCXj) - J x L p(x') dx' - R. (1) 

where is a uniformly distributed set of numbers on the interval 0 < <1. 

Typically in a computer simulation, the values, R^, are supplied by a suitable 
pseudo-random number generator on the interval from 0 to 1. A large sample of 
the variable is then found by solving for x^ In Eq 1 for many random numbers. 

In a Monte Carlo coma simulation the object is to. build up the complete 
distribution of the observed species (atoms or radicals) in space at some 
absolute observation time, t^g. In order to accomplish, this parent 
molecules are emitted from the nucleus during a long time interval of length 
t^, preceeding For a constant production rate (a steady-state model), 

the density distribution function which describes the emission times for 
parent molecules is constant; in other words all emission times are equally 
probable. Equivalently, all emission backup times before the observation 
time, t k , are equally probable. Therefore, as in the original Monte Carlo 
model (Combi and Delsemme 1980a) , individual parent molecules are emitted from 
the nucleus at many backup time intervals, t^, distributed randomly on the 
complete coma buildup time interval, tf. This results in a simple application 



of the Monte Carlo method so that the parent production backup time intervals 
are given simply as tj-t^R^. For an isotropic point source the initial parent 
molecule ejection direction is defined in terms of the standard spherical 
polar and azimuthal angles (0. and as cos 6 ^ - 1 - 2R^ and ^ - 2 jtR^, 

where a different random number is used for each angle. In these simple cases 
the probability density function is constant and the application of the Monte 
Carlo method is trivial. 

Certain processes like those just mentioned can be characterized as an 
even distribution, e.g., for a steady state vaporization rate all production 
times are equally probable. Even for case of dissociation with a constant 
lifetime all population probabilities are equally probable (Combi and Delsemme 
1980a). However, owing to overall poor statistics resulting from the modeling 
of a relatively rare event such as a dissociation, it is desirable to force 
the event to occur in a specified amount of time and penalize the resulting 
particle by only counting a fractionally weighted particle. Similarly, an 
uneven statistical coverage in the output bins, which results when the radial 
space and column density bins are distributed logarithmically in distance from 
the nucleus, can be compensated for by skewing the known probability density 
distribution in time and then by correctly weighting the final particle. In 
our current MCPTM we use both the procedures of forced dissociation and that 
of skewed production and dissociation time distributions to make the model 
much more efficient, and also introduce the inclusion of time -dependence in 
molecular lifetimes. 

In the original steady-state procedure, a decay time interval, tjj, was 
computed using the principle of Monte Carlo, which yields the relation 
tp — -r In (1-R^) where r is the exponential decay lifetime. This was 
computed twice, once for the parent and once for the daughter; a daughter 



trajectory was only computed if the parent decay occurred within the backup 
time interval t^ and the daughter decay did not. A better method is to use 
the concept of a forced dissociation for the parent and a weighting function 
for the daughter decay (see forced first scattering discussion by Cashwell and 
Everett 1959) . To force a parent dissociation within the backup time interval 
t^, the probability distribution is renormalized within that interval, and the 
dissociation time interval is given by 


t 


D 


- r In 


1 




and the particle is assigned a weight (** 1) given by 


( 2 ) 


w D - 


1 


e 


-tiA f 


(3) 


where r is the parent decay lifetime. The weight, w^, of the daughter from 
its own exponential decay is given by 


w . 


-^ t l" t D )/ r d 


<*> 


where is the daughter decay lifetime. 

In the time -dependent model, the variable production rate is treated 
quite simply by weighting (wq) a parent molecule, according to the production 
rate at the emission time (^observation " t i^ relative to the production rate 
at some standard time, e.g., at the observation time or at a comet helio- 
centric distance of 1 AU. Thus a water vaporization curve or a simple power 
law in heliocentric distance can be assumed; for that matter, transient 
activity or outbursts can also be accounted for by simple weighting. 
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The variable lifetime for the parent and daughter decay are somewhat more 
involved. The generally time -dependent lifetimes r - r(t) can be approximated 
as a set of finite element steps 

T “ r 2^2^' .... r n (T n )] 
where (T-^, T 2 , T 3 , . . .T n ) are sub- intervals in time such that 


T 1 + T 2 + T 3 + . . . + T n - ti. 


The partial weight associated with the forced dissociation is then given by 


W D ™ 1 * ex Pt‘( T i/ r i + T 2/ r 2 + T 3/ r 3 + ••• + T n / r n^ 


(5) 


The probability for a dissociation somewhere on this interval 0 < t^ < t^ is 


P(t D ) - 1 - e-r 


( 6 ) 


where 


VV T m-1 V <T 1 + T 2 + • • • + T »-l> 

P — “ + — + . . . + 


r l t 2 


m-1 


m 


( 7 ) 


and 


T j_ + t 2 + • • • + T m ^ t n 5 T, + T 0 + ... + T_ (m < n) 


m-1 D 1 2 



The quantity, p, is just the number of mean free times (like mean free paths) 
contained in time tp. 

Just as in the steady-state case, the random number enters as the 
fraction of parents dissociated until the time tpj therefore, 

R £ - (1 - e'^)/w D (8) 

and solving for p we have 


p In (1 - RjWp) 


( 9 ) 


Finally, tp is determined from the inequality 


Ti T~ T n T, To T_ 

7 1 + T 2 + ... + 7 ^ 




r 1 *-* < p < 7^ + + ,.. + 7“ 

r m-l T 1 r 2 r m 


(10) 


and the equation 


C D “ T 1 + T 2 + 


+ T„ 1 + r_ 
m- 1 m 


Zi 

T i 


To 

+ r + 

r 2 


. + 


r m-l 


(U) 


Once tp is determined, the daughter decay partial weight can be 
calculated in a more straightforward manner by simply summing up the partial 
exponential weights for the remaining intervals from T m up to T n . 


rr 


m 


w^ — exp - 


T 1 r 2 T m 

— + — +...+ - P 

r l T 2 


m 


r m (daughter) 


( 12 ) 


W 


iiitx . j n 

r m+ l( daughter) * * ' r n (daughter) 



The final weight associated with the daughter is found simply by the product 
of the partial weights, WqWpW^, 

The above discussion is appropriate for treating the true probability 
density functions. For the case of modeling the inner coma it is highly 
desireable to skew both the parent molecule production times and all dissoci- 
ation times to small values. In the inner coma model radial distance bins are 
typically distributed logarithmically, and in the application of the inner 
coma model to the question of the computation of photochemical heating it is 
especially important to insure constant statistical certainty in all radial 
bins even those in the innermost coma. In order to understand the conse- 
quences of incorporating a skewed probability density distribution it is best 
to return to the basic probability density distribution function description 
of the Monte Carlo principle Eq.l. In the applications of interest for the 
coma MCPTM (either production time or dissociation time) it is convenient to 
construct a generalized state variable, tj, for the process such that the 
probability density distribution function can be defined as p(x)dx - d 17 and 
the solution for the principle of Monte Carlo (Eq.l), is then simply given by 
o/^i dif-R^, which yields . 

For the case of uniform parent molecule production times the state 

variable is given by qj-tj/t^; for constant dissociation lifetimes the state 

variable is qj— (l-e* t I)/ r )/(l-e" t i/ r ) . When each of these definitions is 

equated to and solved for t^ and tp one finds, as expected, that t^-t^R^ and 

that tp is given by Eq 2 . After trying a number of functional forms we found 

the best way to skew the resulting distributions toward smaller values of t- 

and tp was to require a distribution of the state variable of the form 
r x i 

J p(x)dx - q.-R. , where n is a positive number. It can be shown that a 
o 


17 


normalized skewed probability density distribution function of the state 
variable of the form 


P s (»7)d»? 


1 <h-i> 

h d *» 


( 13 ) 


when integrated from 0 to results in the desired distribution. However, 
since this is not the true flat distribution in the appropriate state vari- 
able, a corresponding fractional weight must be included when the particle is 
counted in the model. Since the true physical distribution is flat then the 
appropriate weight for the skewed distribution of some value of the state 
variable, r / , is found simply by the inverse of the probability density at that 
value. This is 


w 


n>7 



(14) 


The weights for the skewed production time distribution, w^, and the skewed 
dissociation time distribution, w^, are 


w. 


n 


t- d-^) 

C f 


(15) 


and 


w D - n 


1-e 




'■1-e 


-t.A 


(l-f) 


n' 


(16) 


In the final model output bins, each particle is counted as a weight and 
all particles contribute. The total weight is set by the product of several 
partial weights: Wq for the production rate, w-^ for the skewed production time 


distribution (Eq 15), for the forced dissociation (Eq 5 or 16), and for 
the decay of the particle itself (Eq 4 or 12). The value of n in Eqs. 13-15 
can be adjusted to produce the desired even statistical coverage in all output 
bins. A value of 3 has been found to produce reasonable results when a value 
for the total model run time for building the observable coma, tf, is 
typically chosen to be at least 7 times the particle (radical or atom) 
lifetime. 

Most of the discussion to this point has dealt with using the Monte Carlo 
method to compute the weight associated with a particle that results from its 
production and decay mechanisms. There are also Monte Carlo applications 
associated with the trajectory calculation of the particle. In current 
versions of the MCPTM these applications include: (1) speed distributions 
from excess photolysis energy, (2) Maxwell -Boltzmann distributions for the 
convective Maxwellian description of an outflowing coma having a specific 
temperature distribution, (3) multiple collisions of suprathermal particles 
(radicals or atoms) produced by the excess photolysis energy, and (4) random 
ejection directions after dissociation or collision. The method for gener- 
ating random ejection directions was given by Combi and Delsemme (1980a) and 
mentioned again at the beginning of this section. 

Simulating speed distributions (either known or assumed) is a simple and 
direct application of the Monte Carlo method (Eq. 1). In the case of an 
assumed or calculated speed distribution which is imparted on a particle (atom 
or radical) upon photodissociation of its immediate parent (radical or mole- 
cule) there is generally no known analytic form. So a numerical look-up table 
is generated from the speed probability density distribution function, f(v)dv, 
using Eq. 1, appropriately replacing p(x) by f(v), and integrating and solving 
for v^ to form the function, V£-v(R^), in 



(17) 


o J X f(v)dv - R ± . 


An example is illustrated in Figure 6 of Paper 2 for the case of the produc- 
tion of hot H atoms by the photodissociation of water molecules as derived 
from the results of Festou (1981b). In many cases, however, as for CN (Combi 
1980) , and for C 2 (Combi and Delsemme 1986) , we assume an average monoener- 
getic speed. 

Since the outflow of the coma gases is more realistically described by a 
bulk radial speed and some radially dependent temperature we have included a 
convected Maxwellian, and the simulation thereof, as an important mechanism 
in the MCPTM. It is especially important in the case of the photochemical 
heating calculation and the resulting speed distribution of the partially 
thermalized H atoms produced in the photodissociation of water. At some 
point in the coma, the gas flow of a parent molecule can be characterized by 
a temperature, T. The speed distribution, f(v) dv, of molecules of mass, m, 
at temperature T is given by the Maxwell -Boltzmann distribution, 


f(v) 





dv , 


where v m is the most probable speed. 


v m - (2kT/m) 1/2 . 


(18) 


In order to simulate this speed distribution on a microscopic level 
with the particle- trajectory model, we can again use the principle of Monte 
Carlo (Eq. 1). Since the Boltzmann distribution is a normalized probability 


density function, a simulated distribution can be generated from a random 
number sequence, R^ on the interval from 0 to 1 by solving for v^ in the 
equation, 


~Jy 2 / 1 v 2 exp(-v 2 /v m 2 ) dv - Rj . (19) 

7T ' V 0 
m 


The integral in this equation is not analytical, but can be transformed into a 

single parameterized function by making the substitutions: u « v/v m , 

u^ m v i /v m and du - — dv. The equation then becomes 

m 


4 

x 1 / 2 



u 2 e 




du “ Rj * 


( 20 ) 


which has been numerically solved and included as an interpolated table of the 
form, u^ — u(R^) , in the model. 

The final application of the Monte Carlo method to be discussed here is 
that for collisions between the particles (atoms or radicals) and the back- 
ground outflowing coma gas. As discussed in section I, Kitamura, et al. 

(1985) have published the first attempt of performing such a calculation as it 
applies to the thermalization of hot H atoms produced by the photodissociation 
of water molecules. Although we have followed their basic derivation of the 
application of the Monte Carlo method to the collision path, two of their 
assumptions regarding the conditions and density of outflowing coma can be 
improved. One is their assumption of a monoenergetic constant outflow speed 
for the coma. The other is their description of the density in the outflowing 
coma to be simply that for primary water molecules whose population decreases 
both (correctly) as the inverse square of the distance from the center of the 



nucleus and also (but incorrectly) with the exponential scale length for the 
photodestruction of water. 

The first improvement is to consider the coma to be described with a 
variable outflow speed and a variable temperature. Even for rather small 
values for the coma temperature, the thermal speeds of H atoms are generally 
at least comparable and possibly even several times larger than the typical 
values for the outflow speed (~1 km s'^). Kitamura et al. (1985) dealt 
strictly with the effect of the outflowing gas on the final distribution in 
space and velocity of the H atoms. Although this is an important concern, we 
are also concerned with the heating effect of the hot H atoms on the coma. 

They have assumed a form for the description of the outflowing coma and 
performed a Monte Carlo multiple-collision simulation to calculate the H atom 
trajectories. We have found that the MCPTM can be used to calculate the 
photochemical heating which when coupled iteratively to a gas -dynamic model 
yields the outflowing coma description (i.e. the radial dependence in the gas 
temperature and outflow speed) . As is shown in Paper 2 , this coupled approach 
can be used to generate a self-consistent picture of the spatial morphology of 
the H Lyman-o coma and the inner coma conditions. More detailed discussions 
of these issues is presented in section IV of this paper and in Paper 2. 

The second improvement, in addition to providing a better physical 
description of the coma, actually simplifies the collision calculation signif- 
icantly. Kitamura et al. (1985) asstune that the variation with distance from 
the nucleus of the density of molecules from which the modeled particles (H 
atoms in their case) collide is described by the strict density distribution 
of whole neutral water molecules expanding with a constant outflow speed and 
decaying due to photodestruction. Such a description is given by 


( 21 ) 


N H 2 0< r > 


Q 1 -r/(vr) 

4jtv 2 e 
r 


where Q - the water production rate, 
v - the outflow speed, 

r — distance from the center of the nucleus, 
t — the photochemical lifetime of H 2 O, 
and % 2 0 ” t * le number density of water. 


Although this is correct for whole water molecules , this does not describe the 
radial dependence of collisional targets even for the case in question, i.e. a 
pure water coma with a constant outflow speed. As H 2 O molecules are destroyed 
photochemically they produce various fragments, including H, and OH primarily 
but also H 2 , 0 and various ions. All of these fragments also provide targets 
with which the modeled particle can collide. Since, to a very good approxima- 
tion (Johnson 1986) the effective cross sectional area presented by a molecule 
for elastic collisions at low to moderate energies (-2 eV) is given approxi- 
mately by the sum of the individual atomic cross sections, then ignoring the 
photodecay lifetime provides a better approximation to the radial dependence 
gas density as it applies to calculating collision path lengths. Furthermore, 
the absence of the exponential term, as we shall show shortly, provides an 
analytic solution for the probability density integral in the application of 
the Monte Carlo method for the case of hot H atom collisional thermalizatid'n. 

Although we begin with the same formal application of the Monte Carlo 
method as Kitamura et al. (1985), the entire derivation will be presented here 
in the interest of completeness. The first step in this method is to calcu- 
late the collision path length for a single particle having a given location 
in the coma, and a given velocity vector direction. At least the first few 
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general steps also appear in the work of Cashwell and Everett (1959). Since 
the Monte Carlo method is inherently probabilistic, and we wish to derive an 
expression for the collision probability, we begin by considering a beam of no 
particles traveling along a path with a linear coordinate, s, in a stationary 
medium with a number density, N(s) , of scatterers each with a cross section of 
area, a. If the beam of no particles is displaced from s-0 to s-ds then by 
simple attenuation the number of particles scattered out of the beam, dn is 
given by , 


dn - - noN(0)o ds. 


In the case of a medium of uniform density, simple integration of this 
expression yields the exponentially decreasing attenuation law. 

The case of interest here is that for a coma with a highly nonuniform . 
density. For a nonuniform density this equation can be integrated (at least 
formally) to yield 


ln(n/no) - ~a q/ s N(s")ds" . 


( 22 ) 


This can alternatively be written 


-o f S > N(s")ds" 
n _ o J 

n 

o 


(23) 


Returning to the Monte Carlo case of looking probabilistically at a single one 
of the n particles, Eq. 23 yields directly the probability for one particle to 
undergo a collision after traveling a distance, s', in the scattering medium to 


be 
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P(s') - 1 -e 


■o Q f S N ( s " ) ds ” 


(24) 


The probability density function for the collisional scattering process for 
use with the Monte Carlo principle (Eq. 1) then becomes 


, -a / S N(s")ds" 

p(s')ds' - d(<7 Q / s N(s M )ds")e ° 


(25) 


Applying the Monte Carlo principle then yields the relation between a given 
random number and the path length for the collision of a single particle in a 
medium of arbitrary density to be, 


1 - R - e 


-<7 o J s N(s')ds' 


(26) 


which becomes 


a 

o 


J S N(s')ds' 


-In (1-R.) . 


(27) 


In the case of the cometary atmosphere, we shall assume (as discussed 


above) that the coma density is described as 






N(r) 


_Q_ I_ 

4wv 2 ’ 
r 


(28) 


with r measured from the center of the cometary nucleus and, where the 
parameters are defined as in Eq. 21. We shall assume for this discussion that 
between collisions all the modeled particles (atoms or radicals) travel simple 
straight line trajectories. Even when radiation pressure acceleration is 



important, the typical collision path lengths are much shorter than deviations 
of the true path from a straight line. Outside the very inner collision 
region, the radiation pressure accelerations are included in the trajectory 
calculations for relevant species. In situations where the radiation pressure 
is too large to ignore, the curved (parabolic) path can be divided into a set 
of straight line segments, each treated as discussed here. 

The geometry relevant for the calculation of the collision length in a 
cometary atmosphere is illustrated in Figure 2. A particle is located at a 
distance ro from the center of the nucleus and is moving in a direction such 
that its velocity vector makes an angle 0 with the outward radial vector . At 
some arbitrary displacement along this path, I' , the density of molecules in 
the coma (Eq. 28) is given by 

N[r(i')] <r 2 +i» 2 +2r o i , cos«)T 1 . (29) 

However, since the particle is traveling through a radially outflowing coma, 
the path for collisions must be calculated relative to the moving coma gas and 
not in a (quasi-) stationary nucleus -centered frame. The exploded view in 
Figure 2 of the particle located at an arbitrary displacement i' along the 
real path shows the appropriate projection geometry needed to calculate the 
apparent path through the moving gas. A real displacement A i along the real 
path can be separated into radial and tangential components, (Ai)cosfl' and 
(Af)sinfl', respectively. If the particle Is traveling at speed Vg along the 
real path and the coma is flowing radially outward with a speed v, then, 
during the time the particle travels a distance A i along its path, the coma 
gas molecules move a distance Ar - (v/v^)Ai radially outward. The apparent 
path of the particle through the gas (As) is given by the vector sum of the 
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real tangential displacement and the relative radial displacement. This 
yields 


As - Ai[l - 2/Jco s6* + p 2 ] 


1/2 


( 30 ) 


where p is defined as the ratio of the radial outflow speed of the coma to the 
particle speed ( p ■ v / v jg) • 

Since the coma density is a scalar function, it does not depend on the 
relative path between the particle and the coma, but only on the actual 
location in the coma. Therefore N(s) in Eq. 2? may be replaced functionally 
by N[r(s(i'))] « N[r(i')] as given by Eq. 29. A change of variables in Eq. 27 
from ds to di' may then be performed according to Eq. 30. Substituting Eq. 29 
and the differential form of Eq. 30 into Eq. 27 and integrating along the 
apparent path in the rest frame of the nucleus yields the integral equation 


-in(l 
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+ r cos 0 
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+ p 
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di 1 


(31) 
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where 
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+ i' 2 


+ 2r Q i'cos0) 


1/2 


There is no general analytical solution to this integral so that it can be 
solved for £ in the desired form, i - J?(R However, it can be solved in 
certain limits or tabulated through numerical integration. 

One important limit is that for a large particle speed (/9 -♦ 0) , which is 
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important for the case of the collisional thermalization of hot H atoms. In 
fact, the original collision path calculation presented by Kitamura et al. 
(1985) is not the general case derived here (Eq. 31) but is only rigorously 
valid for /9-0, and v - constant since they neglected the relative motions of 
the atoms and the coma gas molecules. By setting to zero, Eq. 31 reduces to 




di ' 


i' 2 +2r f 'cos0+r 2 
o o 


which can be integrated analytically and solved for i for each collision of 
each particle. The result for a single collision path is 


sin0 tan 


(^5q) r o sin * ln < 1 * R i> 


- COS 


'}• 


(32) 


For this case of spherical outflow at a constant speed, the group of 
parameters (Qa/4»rv) which occurs naturally in Eq. 31 is the traditional 
definition of the collision zone radius. It should also be pointed out that 
the argument of the tangent term (the term in []) In Eq 32 is well defined 
physically only for values less than n/2 . As the argument approaches k/2 from 
below, which happens more often for larger values of ro , the collision, path 
length approaches infinity, therefore values out of range imply complete 
escape of the particle. 

For the case of spherical outflow with a variable speed as is implied by 
the results of gas -dynamic models (see section IV) the simple inverse square 
law in density does not exactly hold. In the most general situation Eqs. 28 
and 29 must be modified in order to include a specified speed as a function of 
distance from the center of the nucleus, v(r) . Then the integral in Eq. 31 
would need to be solved numerically. However, the speed is a very slowly 


changing variable in r as compared with any reasonable collision lengths, even 
in the outer coma where the collision paths are large as compared with the 
distance to the nucleus. Therefore, an excellent approximation can be made, 
which is to assume that the collision path is always given by the simple 
inverse square law referenced to the local value of the density. This is 
accomplished by replacing the speed, v, in Eqs. 28 and 32, by v(ro). This 
results in the following expression for the single Monte Carlo collision path 
length £ to be 
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A similar approximation can be made for the more general case, Eq. 31. 

Once the particle is displaced this distance, £, in the coma, it is 
scattered from an outflowing coma gas molecule. To this point we have assumed 
all collisions to be elastic and hard-sphere, where the mass of the molecule 
is the mean molecular mass of the coma gas. The coma at the given location, 
r^, can in general be considered to have a bulk outflow speed, v(r^), and a 
temperature, T(r^) . The outflowing coma molecule is assigned a velocity in 
space, given as the vector sum of the radial outflow speed and a randomized 
thermal component. The random thermal speed component is given by Eqs. 18-20, 
whereas the random direction angles (in the radially outflowing frame of 
reference) are given by cos 0^ -» 1 - 2R^, and - 2 jtR.£, again different 
random numbers are used for each angle. From the three -space velocity vectors 
of the modeled particle (radical or atom) and outflowing coma molecule, the 
straightforward mechanics of an elastic hard-sphere collision is performed. 
First, the center-of-mass of the collision between the given motion of the 
modeled particle and the thermal and radially outflowing molecule is found. 


In an elastic collision, the exit velocity of the incoming modeled particle is 
isotropicized in the center-of-mass frame, given conservation of energy and 
momentum . The isotropization is done by computing two new random ejection 
directions in the center-of-mass frame and then transforming back to the coma 
(lab) frame. Also, the heating of the coma by each collision is calculated in 
radial intervals simply as the mechanical energy transfer per collision, which 
is given by 


AE-|m|Av| 2 (33) 

where Av — change in velocity per collision in the coma (lab) frame, and 
m - mass of the modeled particle (radical or atom). After the new velocity 
vector is computed a new collision path can then be calculated. The modeled 
particle's trajectory is then followed until the the observation time, t Q ^ s . 

In the actual model computer code, modeled particles produced very deep 
in the collision zone (r < 0.01 R c for H atoms, r < 0.1 R c for heavy radicals, 
where R c -K)o/4»rv) , are assumed to be completely thermalized locally. They are 
then transported upward with the outflowing coma gas to the adopted boundary 
distance and re-released with an assumed random thermal plus outflow speed. 
Otherwise the particle could have to be followed for a very large number of 
collisions, which would be a needless waist of time. By choosing an optimal 
boundary the proper physics can be modeled in the most expedient manner. 

Returning briefly to the expression for the collision path (Eq. 32), it 
is evident that there is both a directional ( 0 ) as well as the radial 
dependence to the path length. Similarly to the case of a scattering medium 
with a uniform density, a mean free path can also be derived for the inverse 
square density dependence in the coma. The mean free path is defined such 


that the abundance in a beam of particles travelling through a medium falls to 
1/e of the original abundance. Then for Eq. 32, the mean free path is that 
value of i such that - 1-e"^. In this case the mean free path has both a 
radial and a directional dependence. Figure 3 shows the loci of the mean free 
paths at a number of cometary distances relative to the collision zone radius 
R c . Deep within the traditional collision zone the mean free path for 
collisions is not only small compared with the distance to the nucleus, it is 
also nearly independent of direction. When reaching distances of the order of 
the collision zone radius, the mean free path not only becomes much larger, it 
also becomes highly dependent on direction and favors radial escape. 

In the case of hot H atoms, the large mass ratio between the H atoms and 
the generally heavy coma molecules causes the excess energy transfer per 
collision to be small and thus many collisions are required for the H atoms to 
thermalize completely. By inspection of Figure 3, then, one can see how the 
thermal izat ion efficiency begins to fall from 100% at one-sixth to one-tenth 
of the collision radius, owing to the combination of the large mass ratio, 
large collision paths and favored radial escape. 


III. Particle Trajectories in the Extended Coma 


As discussed in section I of this paper, models of the large extended 
hydrogen cloud have been developed in order to explain the spatial morphology 
of the observed Lyman-a coma (Keller and Thomas 1975, Keller and Meier 
1976). Festou et al. (1979) also applied the vectorial model (Festou 1981a, b) 
to analyze Lyman-a spectral line profiles of comet Bradfield. However, the 
models only addressed the inner coma, and were therefore (appropriately) 
steady- state, one -dimensional models which properly accounted for the 
isotropic ejection of H atoms upon photodissociation of OH and H 2 O. On the 
other hand, the extended cloud models of Keller an co-workers included the 
effects of an isotropic point source of H atoms in which the outflow could be 
described by a Maxwell -Boltzmann speed distribution, solar radiation pressure, 
a variable production rate (typically of the form Q-Qgr’ n ) , and an ionization 
lifetime for the H atoms (of the form r— rgr^) . The effect of the different 
relative heliocentric orbits of the H atoms and the nucleus was accounted for 
by emitting the atoms from the syndyname as projected on the sky plane. 
Heliocentric distance and velocity dependent quantities were evaluated on the 
basis either of conditions at the location of the nucleus or on a column line - 
of- sight average. 

Since the MCPTM method follows individual particles it is then possible 
to model even the three-dimensional, time -dependent extended distributions of 
coma atoms like H, C and 0 without the need of compromising approximations. 

Of course in some instances making approximations and reducing dimensionality 
does make the MCPTM more efficient (as is implicitly assumed in using other 
methods) however it is also possible to verify any approximations explicitly 
by considering a more general case. (This is difficult at best using other 
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methods.) A more detailed discussion of the advantages and numerical details 
appropriate to the H cloud model in particular and the resulting difference 
between the MCPTM and other modeling efforts are given in the companion paper 
(Paper 2) . 

The extended neutral cloud model uses the inner coma MCPTM described in 
the previous section for its source. It is therefore only necessary in this 
section to describe those parts of the extended MCPTM which pertain to the 
large scale geometry, sky plane projection and the orbital trajectory 
calculation. We have thus divided this section of the paper into two sections 
which deal with these areas. Some the numerical details have been placed in a 
separate Appendix in order to facilitate the flow of the derivations in this 
section. 


Solution for the Heliocentric Trajectory of a Particle 

The equation of motion for an atom or molecule in the coma may be readily 
reduced to a modified classical two -body problem of the form 



(34) 


where 


H - GM - S 


(35) 


2 

Here, G is the gravitational constant, M is the mass of the sun, and S/r is 
the acceleration of solar radiation pressure on the gas atom or molecule. The 
simplest case occurs when S is a constant so that (34) reduces to the standard 
two-body equations. For hydrogen atoms, however, S is not a constant, but 
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depends upon the radial velocity of the atom relative to the sun because of 
the doppler shift out of the H Lyman- o absorption feature in the solar 
spectrum 



r* T* 

<r-r) 


(36) 


Because of the dependence in Eq. 36 on the radial velocity of the atom, the 
equation of motion (34) is integrated numerically by adopting a fourth-order 
Runge-Kutta method and specifying appropriate initial conditions for the 
position and velocity of the atom. Equation (34) is solved in the helio- 
centric inertial frame (x,y,z) discussed below (also see Appendix) and 
illustrated in Figure 4. Specification of the initial location (xQ.yg.Zg) 
and velocity (x^.y^.z^) of the atom in the comet plane provides the required 
initial position (x^.y^.z^) and velocity (x^.y^.z^) in the inertial frame 
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that are needed to solve (34). In the above expressions, the location (r c ,f) 
and velocity (V r ,V£-) of the comet at the time of the initial conditions are 
determined by a method to be described below. 

Geometry for the Comet and the Cometarv Atmosphere 

To describe the spatial morphology of gas atoms and molecules in the 
comet atmosphere, a geometric mapping that relates the location of an atom or 
molecule in an inertial sun-centered coordinate frame (x,y,z) to its apparent 
position in the earth sky-plane (M',N # ), centered on the comet, is required. 

In Figure 4, the comet location, comet orbit plane, and comet orbit are 
illustrated relative to the Earth and Sun. Also depicted in Figure 4 are 
the inertial frame (x,y,z) and two comet centered coordinate frames: the 

Sun-oriented (x,y,z) coordinate frame and the historically adopted (see Finson 
and Probstein, 1968) Earth-orientated (L,M,N) coordinate frame. Including 
parallax projection, the Earth sky-plane coordinates (M',N') for an atom or 
molecule at location (L 0 ,M 0 ,N 0 ) in the (L,M,N) coordinate system are defined 
by the parallax projection on the (M,N) plane, and are given by 


M ' - v 1 + tV> 

c o 


(37a) 


N ' “ N o (1 + 

c o 


(37b) 


where p c is the distance of the comet from the Earth. 
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The Sun-centered inertial coordinate frame (x,y,z) is defined with the 
z-axis normal to the comet orbit plane , and with the x-axis directed from the 
Sun through the perihelion point of the comet orbit located a distance q from 
the center of the Sun. The trajectories of atoms or molecules escaping the 
comet may be calculated in this coordinate frame by methods described already. 
The comet location in the (x,y) plane is determined by the angle f between the 
x-axis and the Sun-comet direction, and by the sun-comet distance r c . The 
quantities r Q and f , which depend upon the comet orbit eccentricity e and 
perihelion distance q as well as the difference between the absolute obser- 
vation time t and the perihelion passage time r, are calculated by methods to 
be described later. The (x,y,z) comet centered coordinate frame is defined so 
that the (z-axis) is normal to the comet plane, and the positive x-axis is 
along the Sun-comet direction so that the Sun is located at (-r c ,0,0). A 
point in the inertial frame (x,y,z) is therefore transformed tp the (x,y,z) 
frame by a simple rotation about the z-axis followed by a simple translation 
along the sun-comet direction: 
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(38) 


The (L,M,N) coordinate frame is defined with the positive L-axis along 
the comet -Earth direction, the positive M-axis perpendicular to the L-axis and 
in the plane formed by the Sun-comet and Sun-Earth lines, and the N-axis 
normal to the (L,M) plane as shown in Figure 4. The transformation between 
the (L,M,N) frame and the (x,y,z) frame is given by (see Finson and 
Probestein, 1968) 
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where the unit vectors (l^.m^.nj), (ly,niy,ny) and are defined in 

the Appendix and where (X,Y f Z) are the geocentric equatorial rectangular 

coordinates of the Sun (see the Appendix) . To determine the transformation 

(40), the six orbital elements of the comet (w,0, i,q, e, T) , the obliquity c of 

the Earth, and the absolute time t at which the comet is observed must also be 

specified. For a particular atom or molecule location (x ,y ,z ) in the 

P P P 

(x,y,z) inertial coordinate frame, the desired Earth sky-plane coordinates of 
the particle (M^ ) may therefore be determined by using Eqs. 38 and 39 in 
Eq. 37 and may be symbolically denoted by 


M' 

P 


M p (x p .y p .z p |w,0,i,q,e,r|X,Y,Z,£|t) 


(41a) 


N p “ N p (x p ,y p .z p |w,0,i,q,e,r|X,Y,Z,€|t) 


(41b) 


The equation of motion for a comet in the gravitational field of the Sun 
is readily reduced to the classical two-body problem and exhibits orbits that 
are conic sections with the sun at a focus. The conic section in polar 
coordinates (r,f) (see Figure 4), 




r 


1+e cos f * 


and its exact shape are completely determined by the eccentricity e (e<l 
ellipse; e-1 parabola; e>l hyperbola) and the perihelion distance q (P-q(l+e) 
ellipse; P-2q parabola; P-q(l-e) hyperbola). The position and velocity of the 
comet on its conic section orbit at a given time t depend upon the lapse time 
t-r between that time and the perihelion passage time r. By integrating the 
equations of motion, an expression for t-r involving transcendental functions 
of f can always be determined, but the inversion of the expression to 
determine f as a function of t-r (with the exception of the parabolic orbit) 
is not algebraically tractable. This requires that the expression for t-r be 
inverted numerically by using a simple iterative scheme. For an elliptic 
orbit with an eccentricity in the range 0 < e < 0.95 and for a hyperbolic 
orbit with an eccentricity e > 1.05, standard iterative methods (see Roy, 

1965) are employed to determine the polar coordinate position (r,f) and 
velocity (v r ,Vf) of the comet at a lapse time t-r. For values of the 
eccentricity in the range 0.95 < e < 1.05, the Gauss method (see Benima, 
Chernizak, and Marsden, 1969) is definitely a superior numerical method and 
has been adopted here. 



IV. Photochemical Heating 


Deep within the collision zone mean free paths for intermolecular 

collisions are very small compared with the distance to the nucleus and thus 

o 

are small compared with changes in density (due mainly to the dominant 1/r 
distribution) . In the transition region the mean free path is however highly 
dependent on direction and favors radial escape, as discussed in section 2 
(see also Figure 3). This effect is especially important for the case of 
collisional thermalization of cometary hydrogen which is produced with an 
excess kinetic energy of -2.6 eV during the photodissociation of water. 

Because of the large mass ratio between the heavy parent gas (-18 amu) and H 
atoms (1 amu) only 10% of the excess energy possessed by the H atom is trans- 
ferred to the parent gas per collision. Therefore it takes many collisions for 
an H atom to thermalize completely. For example it takes 7 collisions to lose 
50% its energy, and 22 collisions to lose 90% of its energy. With the direc- 
tional dependence of the mean free path in mind it becomes clear that many H 
atoms produced well inside the collision zone will eventually escape without 
being completely thermalized. Ip (1983) and Crovisier (1984) have recognized 
this fact and have pointed out that the only way of effectively calculating 
this effect is through a multiple collisional Monte Carlo model. The effect 
of the collisional decoupling is important not only for determining the 
photochemical heating rate in the coma but also for determining the outflow 
speed distribution of H atoms from the inner coma which shapes the observed 
morphology of the Lyman-a coma (Paper 2) . 

Both Ip (1983) and Crovisier (1984) have demonstrated that the essential 
physical development (i.e., the radial dependence of gas temperature and 
outflow speed) of a water dominated cometary coma can be reasonably estimated 
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by considering a simple gas -dynamic model of a pure water coma where the only 
heating mechanism is the thermal ization of hot H atoms produced by the photo - 
dissociation of water molecules. In fact, Bockelee-Morvan and Crovisier 
(1987) find that both the gas temperature, as inferred by the near IR obser- 
vation of H 2 O in comet Halley, and the outflow speeds determined from the 
doppler widths of rotational lines of HGN, as measured in radio observations 
of comet Halley agree rather well with the prediction of such a simple 
model. In order to calculate the effect of collisional decoupling of H atoms 
from the coma, a procedure has been developed to iterate between the steady 
state inner coma MCPTM which explicitly calculates the correct heating rate 
and a gas -dynamic model which calculates the outflow speed and temperature. 

o 

Since the approximate 1/r density distribution dominates the collision rate, 
convergence is quite rapid. 

A discussion of the basic procedure as well as the first results of its 
use have already been briefly presented by Combi (1987). We have adopted a 
multi-species single-fluid gas-dynamic model for the physical development of 
the cometary coma. The model is based on the equations of mass, momentum and 
energy conservation and the ideal gas law, which are given as 



( 45 ) 


where p - mass density, Q - gas production rate in molecules per second, 
m - mean molecular mass, v - bulk outflow speed* P - pressure, T - gas 
temperature, 7 — ratio of the specific heats, and S and L are the heating and 
cooling rates. The appropriate quantities, namely m and 7 are species 
averaged (Huebner and Ready 1983) , and the detailed water photochemistry (see 
Paper 2) is followed. It has been demonstrated that the gas-dust interaction 
in the very inner coma may be quite important. Original dusty gas models 
(Marconi and Mendis 1984; Gombosi et al, 1985) found that the main effect of 
dust on the gas dynamics was a mass loading which slowed the gas outflow so 
that inner coma speeds only reached values of -0.65 km/s. More recent work 
suggests the possibility that there is a strong infrared radiation coupling 
between the water -dominated gas and the dust (Marconi and Mendis 1986). This 
results in higher outflow speeds of ~1 km/s and higher gas temperatures of . , 
-200 K at 100 km from the nucleus. 

However, in recent work by Bockelee-Morvan and Crovisier (1987), in which 
they approximated the collisional decoupling of hot H atoms, they have found 
that a simple gas -dynamic model (with a photochemical heating efficiency less 
than unity) can explain the general levels of both the cool temperatures 
inferred from IR spectra of water, as well as the outflow speeds inferred from 
the doppler widths of HCN radio lines observed in comet Halley. . From this 
they conclude that the infrared radiation coupling of gas and dust suggested 
by Marconi and Mendis (1986) is not important and can be neglected. 
Furthermore, the agreement of our modeled Lyman-a isophotes of comet Kohoutek 
(in Paper 2) using only a simple gas -dynamic model also supports this 
contention. 

We begin with the velocity and temperature from a simple pure gas-dynamic 
model. This temperature and velocity field is used as the input to 


characterize the distribution of the outflowing coma gases, as they serve as 
both the extended source for the H atoms and as the collisional scattering 
background for the trajectories of the H atoms: The steady-state inner coma 

MCPTM is then used to calculate the new heating rate explicitly by counting 
the energy deposited in each collision by every H atom. This heating rate as 
a function of radial distance is then used for the term, S in the gas-dynamic 
model (Eq 44) . The resulting temperature and velocity fields from the 
modified gas -dynamic model are then used again in the inner coma MCPTM, and 
the process is iterated until a convergent solution is reached. 

Unlike the original Monte Carlo models (Combi and Delsemme 1980a) the 
coma includes the effects of the variable outflow speed and the temperature on 
the initial space -velocity distribution of the H atoms. The parent dissoci- 
ation time is calculated according to Eq. 2 but now the radial positions of 
the parents are calculated not assuming a constant speed but using the radial 
postion as a function of time, which is generated in the gas-dynamic model. A 
daughter H atom is then emitted from the center-of-mass motion of the parent 
water molecule, which has a radial component and a random thermal component 
from the local gas temperature, with an excess speed determined from the 
excess photodissociation energy spectrum. The details of the water (and OH) 
photochemistry are given in Paper 2. 

The results of a sample calculation of this type are shown in Figure 5. 
The results of the coupled gas-dynamic/MCPTM calculation (lower heavy lines) 
are compared with a simple gas-dynamic model (upper thin lines) for the case 

OQ 

of a 100% water coma at 1 AU from the sun and a production rate of 10 
molecules per second. The pure gas-dynamic results were obtained using the 
standard constant heating rate as in the results of Ip (1983) and Crovisier 
(1984). For both sets of calculations the results are cut off at the 


collision radius since the gas-dynamic model is no longer appropriate at these 
distances . 

The importance of the collisional decoupling of H atoms is most evident 
in the comparison of the heating rates in the two calculations. The H atoms 
begin to decouple fully one order of magnitude inside the traditional 
collision zone, and the rate decreases a little faster than 1/r thereafter. 
This is in contrast to Ip's (1983) estimate of the decoupling which falls 
exponentially since it takes into account the lack of local thermalization of 
locally produced H atoms, but does not account for the collisional heating 
from non- locally produced H atoms. A less severe, exponential form has also 
been suggested by Kitamura (1986) . 

The effect of the hot H atom collisional decoupling is best demonstrated . 
in Figure 3 which shows the directional and radial dependence of the mean free 
path. We thus find outer coma temperatures which are substantially lower 
(30K) than those predicted by simple gas -dynamic models (120K) at the 
collision radius (-3000 km) . The outflow speed is only slightly modified in 
this case, however for comets with a higher gas production (and larger 
collision zone) differences in speed also become substantial (Combi 1987). 

As mentioned the inner coma MCPTM includes the radial dependence for both 
the gas temperature and outflow speed, so the coma gas is considered to have a 
variable convected Maxwellian speed distribution. This is an important 
improvement over the model of Kitamura et al. (1985) which considered the coma 
to have a constant outflow speed and essentially a zero temperature. This 
importance is twofold. First, the heating (collision) rate is coupled in a 
non-linear fashion to the outflow speed and thus is somewhat self regulating 
(through the coupled fluid equations). Second, the exit speed distribution of 
the more thermalized H atoms is a superpostion of the bulk gas outflow speed 


44 


and a (partially) thermalized speed which is substantial for light H atoms 
even at these low temperatures (ten to several hundred Kelvin). This is very 
important for shaping the morphology of the extended hydrogen coma, which is a 
sensitive indicator of the effective outflow speed distribution of H atoms 
exiting the inner collisional coma. An extensive discussion of our early 
attempts to model observed isophobes of the Lyman-a coma of comet Kahoutek 
(Paper 2) demonstrates that the specification of physically realistic condi- 
tions (coma temperature and outflow speed) is required in order to obtain good 
model -data agreement. 

The gas-dynamic model results presented by Bockelee-Morvan and Crovisier 
(1987) are based on a single Monte Carlo simulation for the collisional 
decoupling of H atoms from the heavy outflowing (mainly ^0) gas. They have 
adopted the analytical approximation of Kitamura (1986) to describe the 
results of a single explicit Monte Carlo simulation, and apparently 
extrapolate the analytical form to a variety of cases. They do not explicitly 
include the Monte Carlo calculation in the iterative procedure as we do here 
and as was presented in earlier results (Combi 1987). They, on the other 
hand, explicitly iterate the IR raditive transfer effect in the inner coma 
cooling rate whereas we adopt the approximation of Huebner and Ready (1983) 
for IR trapping. One objection to the analytic form suggested by Kitamura 
(1986) is that the photochemical heating efficiency falls exponentially 
outside some critical radius as did that assumed by Ip (1983) although the 
drop-off is actually not as severe. Through a number of Monte Carlo 
simulations we find, however, that heating efficiency begins to depart from 
complete thermalization much closer to the nucleus than Kitamura' s expression 
but then only falls with an inverse power law in the distance to the nucleus 
having an exponent slightly more negative than -1. 



V . Summary 


In this paper we have presented the mathematical descriptions for the 
methods employed in our general Monte Carlo particle trajectory model (MCPTM). 
These methods allow for the proper treatment to be made for many physical 
processes which are important in both the outer extended atom comae, and the 
inner radical comae of comets. Some of the new methods introduced here have 
been developed in the interest of computational efficiency, while others are 
the inclusion of new physical processes or dimensional generalizations. The 
application of the MCPTM to the spatially extended, time -dependent, and three- 
dimensional Lyman-a coma is the subject of the accompanying paper in this 
volume (Combi and Smyth 1987b, Paper 2). One application which is presented 
in this paper is the proper calculation of the photochemical heating of the 
cometary coma due to the collisions of hot H atoms produced in the photo- 
dissociation of water molecules with the coma gases. In order to study this 
problem the coupled gas-dynamic/MCPTM introduced by Combi (1987) has been 
discussed in detail here. We find that the collisional decoupling of the hot 
H atoms from the heavy parent coma gas well inside the traditionally defined 
collision zone causes the appropriately modeled coma to have lower outflow 
speeds and. substantially lower temperatures. Previous attempts to approximate 
analytical solutions to this decoupling have resulted in the incorrect 
functional form and/or a gross overestimation of the decoupling. 

The methods presented in this paper will be used in future work to 
study the effects on the spatial distributions of observed chemical species 
(radicals and atoms) and to explore the basic physics of the transition zone 
between true fluid-flow to free molecular flow. In our preliminary work, we 
find that the transition zone is large and limits the spatial extent to which 
purely collisional fluid models can be applied. Furthermore, in large, bright 



comets like Halley, typical ground-based observations covered the spatial 
domain of the transition zone itself. Whereas for smaller comets, where 
observations are typically made outside the transition zone, the true initial 
conditions for appropriate free molecular flow (exospheric) models are in fact 
still determined by conditions in the inner and transition regions. 
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Appendix 


The coordinates (x£,yg,Zg) of the Earth in the comet centered coordinate 
frame (x,y,z) were defined in the text by the transformation 


V r c „ 


1- 

m- 

n- 

X 

X 

X 

1- 

m~ 

n- 

y 

y 

y 

1- 

m~ 

n~ 


X 

y 


(A. 1) 


where (X,Y,Z) are the geocentric equatorial rectangular coordinates of the Sun 
or alternatively (-X,-Y,-Z) are the heliocentric rectangular coordinates of 
the Earth. The three rows of the transformation matrix in (A.l) are, respec- 
tively, the three unit vectors of the direction of the x-axis, y-axis and 
z-axis in the geocentric equatorial coordinate system: 


*x “ ^x ,m x ,n x^ 

(A. 2. a) 

iy " 

(A.2.b) 

i z “ 

(A.2.c) 


The geocentric equatorial coordinate system (x,y,z) has its x-axis along the ' 
vernal equinox direction (7) and its z-axis normal to the equator plane. The 
relative orientation of the equator plane, the ecliptic plane, and the comet 
plane are illustrated in Figure A.l. Also depicted in Figure A.l are the 
obliquity of the ecliptic c, the comet plane elements (0,w,i), the sun- 


51 


centered inertial frame (x,y,z), the comet's polar coordinates (r ( 
(x,y) frame, and the comet centered coordinate frame (x,y,z). 

The elements of the unit vectors in (A. 2) follow from Figure 
given by 


*x 


cos (w+f) cos 0 - sin (w+f) sin G cos i 


1- - -sin (w+f) cos G - cos (w+f) sin G cos i 


1~ — sin i sin G 


m£ — [cos (w+f) sin 0 + sin (w+f) cos G cos i] cos e 
- sin (w+f) sin i sin e 




[-sin (w+f) sin 0 + cos (w+f) cos 0 cos i] cos e 
- cos (w+f) sin i sin e 


-sin i cos G cos e - cos i sin t 




[cos (w+f) sin G + sin (w+f) cos G cos i] sin e 
+ sin (w+f) sin i cos c 



[-sin (w+f) sin 0 + cos (w+f) cos G cos i] sin e 
+ cos (w+f) sin i cos e 


n- - -sin i cos G sin e 


+ cos i cos e 


. , f ) in the 
2 and are 

(A. 3. a) 
(A. 3 .b) 
(A.3.c) 

(A. 4. a) 

(A.4.b) 

(A.4.c) 

(A. 5 .a) 

(A. 5 .b) 


(A. 5 . c) 
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Abstract 

The Monte Carlo particle trajectory model (MCPTM) developed in Paper 1 is 
applied to explain the observed morphology of the spatially extended Lyman-a 
comae of comets. The physical processes and assumptions used in the model as 
they relate to the photodissociation of H 2 O and OH and the solar radiation 
pressure acceleration are presented herein. For this first application, the 
rocket and Skylab images of the Lyman-a coma of comet Kohoutek were chosen for 
study. The self-consistent modeling analysis of these data consisted of two 
parts. The first part entailed using a steady-state spherically symmetric 
inner coma MCPTM coupled with a simple gas-dynamic model to. calculate the 
physical development of the coma, i.e. the dependence of coma temperature and 
outflow speed on radial distance to the center of the nucleus, as a function 
of the (time) heliocentric distance of the comet. The inner coma MCPTM was 
used to calculate correctly the photochemical heating of the coma due to the 
partial collisional thermalization of the hot hydrogen atoms produced in the 
photodissociation of water molecules. In the second part of the analysis the 
results from the first part were used in a fully time -dependent and three- 
dimensional extended coma MCPTM which includes the explicit . calculation of 
partial thermalization of the H atoms by multiple collisions with coma 
molecules . 

The same physical model yielded very good matches between the modeled 
Lyman-a isophotes and those observed in both of the two very different images 
of comet Kohoutek. The production rate was varied in time as implied by the 
shape of the visual light curve. All other physical parameters were varied 
only according to their naturally expected heliocentric distance and velocity 
dependencies. The complete physical description of the inner coma provided by 
the coupled gas-dynamic/MCPTM calculation was needed to obtain a good fit to 


the data. The correct inner coma description is important since it provides 
not only the initial conditions for the photodissociated H atoms but also (and 
most importantly) the collisional targets for the H atoms produced in the 
innermost regions -of the coma. Simplistic descriptions for the coma (single 
speed, and perfectly radial molecular motion) do not yield realistic isophote 
contours. The implications of the model results as they apply to other 
comets, species and a variety of conditions are also discussed. 
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I. Introduction 

The time -dependent syndyname models developed by Keller and collaborators 
(Keller and Thomas 1975, Keller and Meier 1976) have been successful in being 
able to reproduce the shape and the radial gradient of the observed Lyman-o 
coma (Meier et al 1976, Meier and Keller 1980). As has been discussed in a 
more general context in the accompanying paper (Combi and Smyth 1987b, here- 
after Paper 1) , the models of Keller and collaborators consider the source of 
hydrogen atoms to be an isotropic point source which follows the syndyname as 
projected on the sky plane in order to approximate the combined effects of 
solar gravity and solar radiation pressure on the atom trajectories. In order 
to reproduce shapes of the two-dimensional sky-plane isophotes as observed in 
resonantly scattered solar Lyman -a emission, they required the use of the 
superposition of two or three Maxwell -Boltzmann distributions for the outflow 
speed of H atoms. Also, the H atoms must have an exponential decay lifetime 
of “2x10^ seconds at 1 AU from the sun. This value for the average lifetime 
can be reasonably well accounted for by the combined decay rates due to charge 
exchange by solar wind protons, photoionization and impact ionization by solar 
wind electrons (with charge exchange being the most dominant) , although the 

I 

actual value can vary considerably with time about the average value (Combi, 

| Stewart and Smyth 1986). On the other hand, it is not so easy to account for 

the required outflow speed distribution. 

It had been demonstrated that, based on circumstantial evidence, the bulk 
of the H atoms and OH radicals observed by their spectral emissions in come- 
tary comae are produced by the photodissociation of water molecules (Keller 
and Lillie 1974) which were originally vaporized from the cometary nucleus. 
Furthermore, the shapes of the visual light curves of comets are consistent 
| with vaporization which is controlled by water ice at the level generally 
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inferred from the H and OH abundances (Delsemme 1973) . The in situ observa- 
tions of comet P/Halley (Balsiger et al 1986, Krankowsky et al. 1986) have 
entirely confirmed this idea. The photodissociation of ^0 molecules by solar 
radiation has most recently been studied in detail by Festou (1981b). It is 
predicted that the photodissociation is highly exothermic, and that most of 
the excess energy is in fact imparted in the form of translational energy to 
the fragments. Nearly 90% of the water dissociated produces H + OH where the 
H atoms have a speed distribution which peaks sharply near 20 and 30 km/s. 

The predissociation of OH by solar photons into the v'— 2 and higher vibra- 
tional levels of the A S state has been studied in detail by Schleicher 
(1983, also see Schleicher and A'Heam 1984) and is believed to be its main 
destruction mechanism. This process is also energetically exothermic and 
should produce H atoms with a speed distribution sharply peaked near 8 km/s. 
More recently, van Dishoeck and Dalgarno (1984) have presented the results of 
ab initio calculations which suggest that absorption of more energetic solar 
photons into higher lying predissociation states result in an additional rate 
which is somewhat smaller but comparable to the rate for the lower lying 
state. The higher energy dissociations are also exothermic and produce faster 
H atoms with speeds of 17 to 26 km/s. • 

The model analysis of two H Lyman-a images of comet Kohoutek taken by 
Meier et al. (1976) on 25 December 1973 and 8 January 1974 when the comet was 
0.18 and 0.43 AU from the sun respectively, implied that the shape and the 
radial gradient of the brightness can be reproduced by an effective outflow 
speed distribution which is formed by the sum of three Maxwell-Boltzmann 
distributions, at 4, 8 and 20 km/s. The existence of speed components at 7 
and 21 km/s had been suggested in older work by Keller and Thomas (1975) on 
comet Bennett (1970 II). Speeds of this order are also generally consistent 


with the doppler widths of conetary hydrogen line profile data (Huppler et al. 
1975, Festou et al 1979). All taken together these results were highly 
suggestive of a major role of H 2 O and OH photodissociation in the production 
of observed cometary H, however many details still needed to be understood. 

In addition, the question of how to best model observations of hydrogen in 
less productive comets in general and all comets at moderate to large helio- 
centric distances remains in question, even if something similar to the 
weighted sums of Maxwell-Boltzmann distributions are reasonable for productive 
comets at small heliocentric distances (like Kohoutek) . 

Simple direct photodissociation should produce a speed distribution 
dominated by sharp peaks and should produce a hydrogen coma having a peculiar 
observable character when the atoms are subjected to solar radiation pressure 
(Keller and Meier 1980) . Furthermore, no direct account can be made for the 
substantial low speed component required for the 25 December image of comet 
Kohoutek and the small low speed component required for the 8 January image. 
Meier et al. (1976) suggested that at least the existence of a low speed 
component could be explained by the thermalization of those .H atoms produced 
well inside the innermost cometary com a. In fact, Kitamura et al (1985) 
demonstrated with a simple steady-state multiple -collision Monte Carlo model 
that, for a productive comet at a small heliocentric distance, a measurable 
thermalized component should result. 

In this paper we present the results of our reanalysis of the two Lyman-a 
images of comet Kohoutek taken by Meier et al. (1976) with our three-dimen- 
sional time -dependent Monte Carlo particle trajectory model (MCPTM). The 
MCPTM includes for the first time a truly physically realisitic description of 
the detailed production mechanisms and trajectories of H atoms produced by the 
photodissociation of H 2 0 and OH. The MCPTM takes into account in a self- 



consistent fashion the combined effects of the multiple collisions of H atoms 
with the outflowing coma gas and the trajectories of the atoms in the presence 
of solar radiation pressure and solar gravity. The results presented in this 
paper represent only the most recent step in an evolutionary process that we 
necessarily followed in order to develop a physically realistic model which 
could reproduce the observed 2-D isophotes of the Lyman-a coma. This process 
should at least be briefly described here in order to provide a clear frame- 
work for understanding the current level of complexity in our modeling 
efforts. 

As discussed above, the past attempts to model the observed distributions 
of cometary hydrogen have fallen into two classes: steady- state inner coma 
models (e.g. Festou et al. 1979, Kitamura et al. 1985), and time -dependent and 
(quasi-) three-dimensional extended coma models (Keller and Meier 1976). Our 
broadly based goal has been to use the power and flexibility provided in the 
Monte Carlo method, in order to combine a physically realisitic inner coma 
model, generalized to include time dependence, with a true three-dimensional 
extended coma model in order to establish a clear physical basis for the 
spatial morphology of the observed Lyman-a coma. 

The first step in this process was to merge the simple Monte Carlo model 
developed by Combi and Delsemme (1980a) , which could describe the photo- 
chemical production of hydrogen from 1^0 and OH with the correct isotropic 
ejection of daughter fragments upon dissociation, with the three-dimensional 
time -dependent trajectory calculation and comet-atom projection geometry, 
appropriate for hydrogen atoms traveling under the combined influence of solar 
gravity and radiation pressure. The extended coma portion of the model alone 
represented a technical improvement over the syndyname model of Keller and 
Meier (1976) in that the physical processes such as lifetimes and emission 



rates did not have to be given average values for a column or at the position 
of the nucleus but could be included explicitly for each atom at each location 
in space and for the entire time history of the development of the coma. 

A collision routine based on a simple single-speed perfect radial outflow 
for the coma gas, similar to that of Kitamura et al. (1985), was then incor- 
porated into the inner coma portion of the model in order to account for the 
partial thermalization of hydrogen atoms by multiple collisions with the 
outflowing coma gas. The inner coma calculation of course had to be gener- 
alized to account explicitly for the time -dependent lifetimes (scale lengths) 
and production rate (see section II in Paper 1), as well as the Delsemme 
(1982) law for the heliocentric distance dependence of the parent molecule 
outflow speed. Once leaving the inner coma the trajectories of the H atoms in 
the extended coma were then calculated explicitly in three dimensions. 

Although the steady- state inner coma Monte Carlo models of Kitamura et 
al. (1985) seemed to account roughly for the necessary fraction of low speed 
hydrogen atoms they made no attempt to reproduce particular observed isophote 
contours and only presented very crude speed distribution functions. In our 
first attempts to reproduce the real observations of Meier et al. (1976) we 
found that the simplisitic description of the coma consisting of a single- 
speed perfectly radial outflow for the coma molecules necessarily produced 
thermalized hydrogen atoms having the same perfect radial outflow speed as the 
coma (l-2km/s using the Delsemme velocity law) . These extremely low speed H 
atoms produced an extremely narrow extended tail not seen In the observations. 
The speed distribution functions of Kitamura et al. were only computed with 
resolution of 5 km/s so the dominance of the extremely low speeds was not 
apparent. A contributing factor which had been overlooked was that even for 
fairly modest coma temperatures (<100 K) the thermal speed component for H 



atoms is at least comparable to the radial outflow speed. For the higher 
temperatures expected for comet Kohoutek at 0.18 and 0.43 AU (for the two 
observed images) the thermal speed components should be much larger. It 
became apparent that some accounting for not only the outflow speed of the 
coma but also for the coma temperature must be included in a collisional model 
for the production of H atoms. 

The first attempt at including a coma temperature was to assume that for 
each of the two images of the comet a single effective time -independent 
temperature could be used for the coma gas. At each collision the target coma 
molecule was given both its radial speed plus a randomized thermal component 
(randomized in both direction and across the Maxwell -Boltzmann distribution as' 
discussed in section II of Paper 1). This simple picture was a tremendous 
improvement in the ability to reproduce realisitic looking isophotes for both 
images. It was therefore clear that the physical state of the inner coma 
(temperature and outflow speed) was an important shaping mechanism for the 
Lyman-a coma. In fact, it was apparent that for comets where collisions 
between H atoms and the coma are important (productive comets at small 
heliocentric distances) that the shape of the hydrogen coma was in fact 
diagnostic of conditions in the inner coma. 

Ambiguities unfortunately remained since it was clear that the coma tem- 
perature varied both with time (heliocentric distance) and with distance to 
the nucleus and since reasonble reproduction of the observed Lyman-a isophotes 
could be made by a range of combinations of the outflow speed and temperature 
for the coma. Looking to the published results of the hydrodynamic models for 
the inner coma (Marconi and Mendis 1982, 1983, 1984, Huebner and Keady 1983, 

Ip 1983, Crovisier 1984) proved to be of little help. Only very few cases had 
been studied. These were usually for productive comets at 1 AU from the sun 
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and no systematic study of the heliocentric distance variations had been pub- 

I 

lished. Furthermore, because of its collisional effect on hydrogen atoms, the 
region of the coma, which is of the utmost importance to characterize, is the 
transition region from collision-dominated flow to free collisionless flow. 
This is the region where the hydrodynamic models are no longer applicable. 

This point had been first emphasized by Ip (1983) who pointed out that 
the main heating mechanism for a water dominated coma is the collisional 
thermal izat ion of the hot H atoms produced in the main (90%) photodissociation 
branch of water. Owing to their small mass many collisions between H atoms 
and the heavier coma gas molecules (H 2 O, CO, OH, 0, etc.) are needed in order 
to thermalize completely. Thus the so-called photochemical heating of the 
coma will become inefficient even well within the traditionally defined 

collision region. Although Ip (1983) suggested an analytic approximation to 

l 

this collisional decoupling, both Ip and later Crovisier (1984) suggested that 
the best way to treat this problem is with an explicit Monte Carlo collisional 
model. Since a collisional routine for the inner coma model had already been 
developed it was possible not only to calculate the effect of the collisions 

on the distribution of a daughter species but also to calculate explicitly the 

' 

energy transfer at each collision to the coma. This energy transfer for the 
many collisions modeled is in fact the photochemical heating calculated at the 
kinetic level. 

It became evident that there were important Interrelations between (1) 
collisional thermalization and decoupling of hot H atoms .(2) photochemical 
heating and physical development of the coma and (3) the effective (and time- 
dependendent) speed distribution of H atoms leaving the inner coma. An 
j iterative scheme was therefore developed that coupled the collisional steady 

state inner coma MCPTM which calculated correctly the photochemical heating 
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efficiency to a simple gas-dynamic model which calculated the outflow speed 
and temperature of the coma. Preliminary model results for comet Halley using 
the coupled gas-dynamic/MCPTM were presented by Combi (1987), however, a more 
detailed description as well as a more general application is presented in 
section IV of Paper 1. 

The results of the last iteration of the general study of the spatial 
morphology of the hydrogen coma, prior to the results contained in this paper, 
were presented in preliminary form by Combi and Smyth (1987a) as the first 
self-consistent calculation which included a time -dependent description of the 
physical state of the coma. The gas -dynamic/ MCPTM procedure was used to 
calculate the coma temperature and outflow speed for a number of heliocentric 
distances of comet Kohoutek. These results were used to infer a set of single 
effective coma outflow speeds and temperatures in the collision transition 
region where the last several collisions occur between the H atoms and the 
outflowing coma. The fact that the coma temperature and outflow speed depend 
not only on time but also on the distance to the center of the nucleus was not 
included. In this attempt very good reproduction for the image recorded on 
January 8, 1974 (r - 0.43 AU) was obtained, however the fit of the image of 
December 25, 1973 (r — 0.18 AU) was not quite as good. 

As will be presented in section IV of this paper, the coma temperature in 
the region where the hydrogen atoms undergo their last few collisions, (about 
10^ km in cometocentric radius) when the comet is near perihelion, is charac- 
terized by a fairly steep drop with increasing distance. Owing to the prefer- 
ence for radial escape, even from fairly deep within the traditional collision 
zone (see Figure 3 in Paper 1), compounded by the large number of collisions 
required to alter the speed of the light H atoms, we found that using a single 
average temperature in this region does not produce the correct effective 


speed distribution for the partially thermalized H atoms exiting the inner 
coma. Furthermore, the outflow speed of the coma, which influences the number 
of collisions through its effect on the gas density (through the continuity 
equation) and also through the dissociation scale lengths (i.e. how far out 
into the coma the H atoms are produced) , also varies with distance from the 
nucleus. It became clear that a completely time-dependent description of the 
whole radial dependence of both the outflow speed and temperature of the coma 
would be required to yield a true self-consistent calculation. 

The remaining sections of this paper deal with the specific application 
of the MCPTM tools developed in Paper 1 to the production of cometary hydrogen 
and to the interrelated aspects of the physics of the inner coma and the 2-D 
spatial morphology of the extended Lyman-a coma . Section II of this paper 
describes the detailed physical processes relevant to the production of 
cometary hydrogen as included in the MCPTM. This discussion deals with the 
photochemical chain for the production of cometary hydrogen atoms, the solar 
Lyman-a emission, collisional cross sections, and the variation of production 
rate with heliocentric distance. Section III describes the heliocentric 
distance dependent inner coma description for comet Kohoutek as modeled with 
the coupled gas-dynamic/MCPTM procedure discussed in Paper 1. In section IV 
the application of the MCPTM to the two-dimensional sky-plane images of Meier 
et al. (1976) are presented and are compared with the earlier modeling 
efforts . Also discussed in section IV are the implications of the success in 
modeling these data on the appropriate models to be used for analyzing data 
from other comets and other cometary conditions. Finally, section V presents 
a brief summary of results contained in this paper as well as a discussion of 
further improvements in and other applications of the current MCPTM to be 


undertaken in the future. 



II. Physical Processes Included in the MCPTM 
Owing to the inherent flexibilty afforded by the many-particle Monte 
Carlo simulation, it is possible to include explicitly the detailed space, 
time and velocity dependencies of the physical processes relevant to the pro- 
duction, kinematics and decay of hydrogen atoms in the cometary environment. 
Separate discussions are presented in this section of the paper regarding the 
photochemistry of 1^0 and OH, the solar Lyman-a line profile and absolute 
flux, the composition and the variation of gas production rate of comet 
Kohoutek, the collisions between atoms and molecules, and the decay lifetime 
of hydrogen. 


Photochemistry of H 2 O and OH 

The photochemistry and expected photochemical kinematics associated with 
the production of hydrogen in comets has been studied in detail by several 
other investigators. The purpose of this discussion is simply to document . 
these details as adopted in the MCPTM analyses presented in this paper. 

Festou (1981b), as mentioned earlier, has studied in detail the various 
branches of the photodissociation of H 2 O by the solar ultraviolet radiation 
field. In addition to assessing the rates of the various branches he has 
calculated the excess energy left after the exothermic photodissociation which 
is shared in the form of kinetic energy between the fragments. According to 
the results of Slanger (1982) a minor revision to the branching ratio between 
the OH + H and the 0(^D) + H 2 for the solar Lyman-a contribution is 
required. Otherwise we have adopted the results of Festou (1981b) to describe 
both the rates as well as the speed distributions for the fragments. These 
results are presented in Table 1. The largest component to the direct source 
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of H atoms at 20 km/s is actually the narrow Gaussian distribution at 19.6 
km/s as determined by Festou, and is included as such in the model. 

The main photodecay process for OH radicals is actually a predissociation 

O . 

from the ground state to the v'-2 and higher vibrational levels of the AS 
state. The transition to the v'-l vibrational level results in the familiar 
observed fluorescent transition at 3090 A. The predissociation process was 
examined in great detail by Schleicher (1983, also Schleicher and A'Heam 
1984) , who calculated not only the dissociation rate but also the heliocentric 
velocity dependence of the rate. This velocity dependence results from the 
Swings effect, since the solar spectrum is dominated by many absorption lines 
(most due to solar OH) whose positions vary relative to the cometary OH 
rotational lines with different doppler shifts. An additional predissociation 
process has been suggested in the results of van Dishoeck and Dalgarno (1984) 
whose ab initio quantum mechanical calculations have identified the existence 
of higher lying states which can be populated by the absorption of solar 
•Lyman-a. They have combined their results for these newly identified 
transitions with the earlier work of Schleicher to present a complete picture 
of the photodissociation rate of cometary OH as well as the energetics of the 
H and 0 atoms produced. 

The absorption of the long wavelength solar UV photons into the A^E + 
state should produce H atoms having excess speeds primarily of 8 and 11 km/s, 
however the higher lying states should produce much faster H atoms with speeds 
of 17-26 km/s. In addition, the results of van Dishoeck and Dalgarno (1984) 
also predict the production of a large number of 0 atoms in the ^D state which 
should produce a much larger spatially extended (although faint) emission of 
the familiar 6300 A forbidden line that has typically been observed as 
concentrated in the inner coma and attributed to one of the primary 
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photodissociation branches of H 2 O (see Table 1) . In recent work by Roesler et 
al. (1987), the emission of the 6300 A line at nearly 10** km from the nucleus 
of comet P/Halley at a level consistent with all of the projected branching 
ratios may be a confirmation of the ab intio calculations of van Dishoeck and 
Dalgamo. Table 2 shows the photodissociation rates and excess speed distri- 
butions adopted in the MCPTM for the OH photodissociation. The expected 
heliocentric velocity dependence of Schleicher (1983) for the A^2 + absorption 
has been explicitly included. 

The speed distributions for the H atoms and OH radicals from all 
photodissociation branches in Tables 1 and 2 are spanned in the MCPTM as 
discussed in section II of Paper 1. The probability distribution function of 
the excess speed was constructed for each photodissociation process , so that 
given* a random number on the interval from 0 to 1 a proper speed for an OH 
radical or an H atoms is chosen according to Eq. 17 in Paper 1. The photo-, 
dissociation rates contribute to the decay lifetimes for H 2 0 and OH and are 
included in the model with their full time (heliocentric distance) dependences 
according to the procedure given by Eq. 5-12 in Paper 1. . The heliocentric 
velocity dependence for the OH lifetime is evaluated according to. the time 
dependence of the heliocentric velocity of the nucleus, since the changes in 
the lifetime are small compared with the dispersion of the speeds of the OH 
radicals in the nucleus (1-2 km/s) . This effect could be treated with the 
MCPTM, but is a needless complication at this point. 


The Solar Lyman-a Flux and Line Profile 
The solar Lyman-a radiative flux serves both as the primary source of the 
radiation observed from the hydrogen coma as it is resonantly scattered from 
the atoms, and (through this scattering) as the source of the radiation 



pressure acceleration felt by the H atoms as they move in their orbits away 
from the comet around the sun. In fact, depending on the level of the 
variable solar flux and on the doppler shift of a cometary H atom relative to 
the solar line profile, the radiation pressure acceleration on H atoms ranges 
from somewhat below to somewhat above the solar gravity acceleration. The 
proper estimation of both the solar Lyman-a line profile shape and also the 
total integrated line flux is of the utmost importance in modeling both the 
production rate of hydrogen (through the photon emission rate or g- factor) as 
well as the speed distribution of atoms exiting the inner coma (through 
radiation pressure acceleration which causes the observed spatial distortion 
in the shape of the isophote contours) . 

The estimate of the Lyman-a flux at line center adopted in the early 
cometary hydrogen work (Keller and Thomas 1975, Meier et al. 1976) of 3.7x10^ 
photons cm’ ^ s'^ for the solar minimum conditions appropriate in the 
period of 1973-1974 is about 40% higher than current estimates of the 
appropriate value (Lean and Skumanich 1983) . There has been some controversy 
over the absolute calibration of different instruments during the last ten 
years although the lower values seem more reasonable. We also now know that 
the actual solar flux in Lyman-a can vary by ±10% over a single 27-day solar 
rotation period. Although it is now possible to assess the day-to-day changes 
in the value by appropriate time shifts (solar rotations) of Solar Mesospheric 
Explorer measurements there are no such data available for the period covered 
in the Kohoutek observations. Therefore, we are required to adopt the new 
average estimated value (2.75 x 10^ photon cm'^ s’ 1 A’ 1 for the line center 
flux and 2.5 x 10^ photons cm’^ s’* for the integrated flux) and to keep in 
mind that variations of the order of 10% may effect the overall production 
rate as well as the radiation pressure acceleration. 
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The solar Lyman-a emission is concentrated in a line whose FWHM is of 
order 1 A and whose line center is self-reversed by about 35% from the maximum 
flux. A full solar disc Lyman-a profile was constructed by Lemaire et al. 
(1978) from an observed high resolution disc -centered and limb spectra. This 
spectrum is reproduced in Figure 1 with the original geocoronal absorption 
line filled- in to approximate the line as seen by a hydrogen atom in a 
cometary coma. In the MCPTM the resonance scattering rate for a hydrogen atom 
is calculated given the doppler shift for the appropriate heliocentric 
velocity and where the integrated line flux, as shown in Figure 1 is scaled to 
the desired value. For example, in the model results published for comet 
P/Giacobini-Zinner (Combi, Stewart and Smyth 1986) the integrated solar line 
flux determined from the measured SME data was used to scale the entire line . 
profile. For the results presented in this paper the integrated flux of the 
Lyman-a line constructed by Lemaire et al. (1978) from observation has been 
scaled down linearly to the adopted value . He're too there has been some 
disagreement in the literature as to whether the integrated flux scales 
linearly with the flux at line center or with some power (Vldal-Madjar 
1977). However, long term observations of interplanetary hydrogen do support 
the linear scaling (Shemansky 1986). The actual calculation of the photon 
emission rate and radiation pressure acceleration from a given solar flux is a 
well-known relation (see e.g. Keller and Meier 1976) and will not be discussed 
here . 


The Gas Production Rate Variation and Composition of Comet Kohoutek 
For the purposes of the modeling analysis presented in this paper the 
adoption of both the variation in time of the water production rate as well 
as a crude estimate of the basic gas composition is necessary. As has been 
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discussed in the past by many investigators, time dependences in the pro- 
duction rate, decay lifetimes and the solar Lyman-o flux likely effect the 
observed morphology of the hydrogen coma in a convolved manner. Without 
independent measures of the day-to-day variations in the solar UV flux, and 
the solar wind flux it is not theoretically possible in a modeling context to 
separate the intertwined effects of variable lifetime and variable production 
rate, although estimates of the longer term averages of these variables can 
certainly be obtained. The composition is important for the coupled gas- 
dynamic/MCPTM calculation of the coma dynamics since the mean molecular mass 
of the coma gas must be estimated not only for the H 2 O which eventually 
produces the H coma but also for the likely presence of CO and CO 2 in non- 
negligible amounts . The contribution to the total density is also needed for 
the collision calculation for the H atoms in the extended coma MCPTM, since 
the addition of 20% of species other than water will decrease all collision 
paths by about 20%, increasing the heating efficiency over a pure water coma. 

In the hydrogen models of Keller and Meier (1976), one of the adjustable 
parameters in the data fitting procedure is the form of the gas production 
rate time variation law. Meier et al. (1976) assumed that the gas production 
rate, Q, varied in time with the heliocentric distance, r, according to the 
commonly assumed power law, Q-Qq 1 "* w * iere Qo production rate at 1 AU 

and n is an adjustable parameter. For each observation (which included the 
two full images we address here, in addition to a set of poor quality images 
that were reduced to upsun and downsun radial profiles) they computed a value 
for Q and n. The upsun and downsun radial profiles are not good indicators of 
the correct outflow speed distribution in the model of Keller and Meier. As a 
result they published a set of H production rates which provided, in effect, 
the heliocentric distance dependence of the production rate. However, owing 
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to large calibration uncertainties in the larger set of poor quality images, 
as well as some of the parameter ambiguities, it would be difficult to draw 
firm quantitative conclusions regarding the variation of the gas production 
rate. Furthermore, apparent variations in the H production rate could have 
been due rather to real day-to-day variations in the solar Lyman-a flux. 

Another approach, which we have adopted here, is to examine the variation 
in gas production rate as inferred from the visual light curve of comet 
Kohoutek. Delsemme (1976) published an analysis, in which he fitted a 
theoretical vaporization curve to a compilation of many visual magnitude 
estimates. For the purpose of this study, the reduced magnitudes of 
Delsemme 's compiled pre- and post? perihelion light curves have been converted 
to the relative time variation of the gas production rate. Figure 2 shows the 
smoothed results of our adopted production rate variation normalized to the 
production rate at 1 AU pre -perihelion, Qq. During the build-up time for the 
hydrogen coma observed in the two Lyman-a images (.43 AU pre -perihelion to. 43 
AU post-perihelion) , the variation in production rate has a power law slope in 
the range of 2 to 3 with a noticeable post-perihelion dimming. The broad 
general trend is at least somewhat similar to the variation inferred by Meier 
et al. (1976) from bulk of the poorer quality upsun and downsun Lyman-a 
profiles. 

This assumed gas production rate variation had been reasonably well 
reproduced by Delsemme' s (1976) vaporization theory which assumed a volatile 
composition of 85% water and 15% of a more volatile substance. This is not 
unlike the composition which can be inferred for the in situ and remote 
observations of comet P/Halley (Krankowsky et al 1986, Balsiger et al. 1986, 
Stewart 1986, Feldman et al. 1986). With all of this in mind we have adopted 
a gas composition for comet Kohoutek which is 80% 1^0, 16.5% CO and 3.5% CO 2 . 


The main effects of the assumed composition are in the mean molecular mass for 
the coupled gas-dynamic/MCPTM and in the coma gas collision cross section in 
both the former and the extended coma MCPTM. 

Hydrogen Lifetimes and Collisions 

The last two physical processes of importance in shaping the distribution 
of hydrogen in the extended coma are the decay lifetime of hydrogen atoms in 
the coma and the mechanical description of collisions between atoms and 
molecules. As discussed in section I of this paper hydrogen atoms in the 
interplanetary space around the cometary nucleus may be effectively removed 
from the coma by three main processes : 80% by charge exchange by solar wind 
protons, 15% by photoionization by solar UV photons, and 5% by electron impact 
ionization by solar wind electrons. Although the canonical value of 2 x 10 6 
seconds for the total lifetime of hydrogen can be derived from mean solar 
extreme UV fluxes, and solar wind density and bulk flow velocity, these 
conditions are known to vary widely on time scales short comparable to the 
lifetime itself. The lifetime of hydrogen atoms has been calculated, using 
data recorded by the ICE satellite during the 42 -day period prior to its 
encounter with the tail of comet P/Giacobini-Zinner, by Combi, Stewart and 
Smyth (1986) for the purpose of the MCPTM analysis of Pioneer Venus observa- 
tions of the Lyman-a coma. It was found that although the 42 -day average was 
nearly identically equal to the canonical value, the running 12 -hour averages 
of the lifetime varied over a range from as small as 2 x 10** up to 3.5 x 10^ 
seconds . 

For the purpose of analyzing the comet Kohoutek data (or that of most 
other comets for that fact) , it is impossible to reconstruct such a detailed 
picture of solar wind conditions needed to calculate the hydrogen atom 



lifetime. Therefore as in the earlier studies we are forced to adopt the 
constant average value. What is important to note however is that, when 
assuming a single constant value for the lifetime in a model (scaled 
appropriately with heliocentric distance) , real variations in lifetime are not 
easily distinguished from variations in production rate. The real lifetime 
variations can be large, and the departures from the average last long enough 
to cause easily detectable effects on the observed abundance and distribution 
of hydrogen atoms. 

As discussed in section II of Paper 1, both the inner coma and extended 
coma MCPTM explicitly calculate the effects of many collisions for the 
daughter H atoms and also OH radicals with the background putf lowing coma 
gas. Each collision has been assumed to be a hard-sphere elastic collision. 
Although this is a fairly simplistic description, it is reasonably appropriate 
for these low to moderate energy collisions between atoms and molecules 
(Johnson 1982). The. collision cross sections between atoms or molecules may 
be well approximated by the sums of the individual atom cross sections as 
determined from the individual atomic radii, which are 0.6 A for 0 atoms and 
0.7 A for H atoms (Allen 1973). These estimates yield total collision cross 
section of 1.89 x 10"^ cm'^ for H onto H 2 O, and 3.24 x 10 cn"^ for OH onto 
H 2 O. These are comparable to the typical values quoted in the literature in 
discussions of cometary collision zones (Festou 1981a, Ip 1983, Crovisier 
1984) . 

For these first results presented in this paper and Paper 1, we have 
assumed all collisions occur between the modeled daughter species (H for 
example) and a molecule having the mean molecular mass resulting from the 
assumed chemical composition. This is a very good approximation if in the 
region where many collisions occur most of the parent molecules have not yet 


dissociated, which is in fact true except in extreme cases of large gas 
production rate and small heliocentric, as was illustrated in Figure 1 of 
Paper 1. It is only when comet Kohoutek was near perihelion that we approach 
such a limit. However, as the results in the next section of the paper show, 
since collisions were dominant out to large distances from the nucleus the 
pure gas -dynamic description is applicable for Kohoutek near perihelion and 
all species are collisionally coupled to one another anyway. Therefore, our 
description of the outflow speed and gas temperature remains valid, and since 
the collision cross sections are calculated as the sums of the individual 
atomic cross sections the collision path lengths are calculated correctly in 
the MCPTM. 

The only uncertainty which remains then is that for the time near 
perihelion the H atoms may have been thermalized a bit too rapidly by not 
accounting for collisions between an appropriately partitioned set of 
species. The inclusion of this effect would necessitate the modeling of the 
detailed kinetic theory for individual species (H 2 O, OH, H, 0 etc.) and would 
entail a tremendous complication of the model far beyond the scope of these 
papers. On the other hand, the general question of the collisional decoupling 
of individual species from each other and of particles in one species from 
themselves represents a fundamentally important problem which the general 
MCPTM method as derived in Paper 1 can ideally address at the kinetic theory 
level. This question will be left for future work. 


III. The Inner Coma of Comet Kohoutek 


Using the iterative approach developed and discussed in section IV of 
Paper 1, the coupled gas-dynamic/MCPTM was used to calculate the physical 
development (outflow speed and gas temperature) of comet Kohoutek at various 
times previous to the two Lyman-a observations which was the goal of the 
modeling analysis. This type of analysis would provide the necessary time 
history important for calculating the correct trajectories for H atoms in the 
coma of the comet. The gas-dynamic/MCPTM uses the steadv-state inner-coma 
MCPTM (see section IV of Paper 1) to calculate explicitly the photochemical 
heating rate of the coma by the most dominant process (by far) which is the 
collisional thermalization of hot H atoms produced by the primary photo- 
dissociation branch of water. A steady-state model for the inner coma can be 
used since the dynamical time scale important for building up the collision 
and transition zone region of the coma (about 10"* seconds) is small as com- 
pared with the time scale involved in changes in cometary conditions due to 
changes in heliocentric distance . These heating rates are then used in a 
single -fluid multi -species gas -dynamic model which calculates the outflow 
speed and temperature of the coma. The .new outflow speed and temperature 
profiles are then used in the MCPTM to calculate new heating rates, and the 
procedure is iterated until a convergent solution is reached. Since the 
collisional inner coma MCPTM is the more computationally intensive of the two 
procedures, the first iterations are made with a low statistics MCPTM (fewer 
total particles) . As the iteration proceeds the number of molecules is 
increased. 

As discussed in section I of this paper the final versions of the Lyman -a 
images as modeled here are the result of an evolutionary process which 
involved less complicated (and less self-consistent) modeling attempts. 


The gas production rate variation was adopted from the visual light curve of 
comet Kohoutek (see previous section and Delsemme 1976) . Both the details of 
the collisional thermalization efficiency in the coupled gas-dynamic/MCPTM and 
the final outflow speed distribution of H atoms exiting the inner coma in the 
extended coma MCPTM obviously depend on the total gas production rate. The 
larger the total gas production rate is, the more efficient the collisional 
thermalization will be and the more H atoms will be partially thermalized and 
slowed from their initial photodissociation speeds. However, from previous 
reasonably successful reproductions of the 2 -D isophote shapes of the January 
8, 1973 image of the comet (r - 0.43 AU) , the appropriate value of the 
normalization factor for the time dependence of the production rate, Qq, in 
Figure 2 was fairly well specified. 

Therefore, for the results presented here the appropriate time dependence 
of the gas production rate could be adopted from this value. Table 3 shows 
the adopted parameters for five gas-dynamic/MCPTM analyses of the physical 
development of the inner coma of comet Kohoutek which span a 16 day period 
prior to the image recorded on January 8, 1974, and cover the required backup 
time to generate the whole Lyman-a coma for both observations. We chose the 
perihelion time for one of the times as well as two pairs of heliocentric 
distances on either side of perihelion, one set by the heliocentric distance 
at the later of the two images (0.43 AU) and one intermediate (0.25 AU) . It 
is not possible to use redundant results for the same heliocentric distances 
on either side of perihelion because the production rate curve is higher 
before perihelion than after and the photochemical heating efficiency depends 
upon the absolute production rate. Furthermore, having results for the same 
heliocentric distance but different production rates provides an interesting 
comparison of the theoretical results themselves. 



The total gas production rate was set for a self-consistent production of 
the observed Lyman-a emission according to the chemical composition adopted 
for comet- Kohoutek and the branching ratios of the hydrogen atoms produced in 
the photodissociation of H 2 O and OH. The total gas production rate is then 
about 36% larger than the value usually taken to be one half of the hydrogen 
atom production rate, which implicitly assumes a 100% water composition and 
100% of the water produces H 4 - OH. Owing to their long photodissociation 
times, CO and CO 2 cannot contribute substantially to the heating rate within 
the collision zone (Ip 1983). The details regarding the functional form of 
the parameters (the IR cooling rates, the radiation trapping) used in the gas- 
dynamic portion of the calculation have been discussed elsewhere (Combi 1987, 
and Paper 1) and will not be reproduced here. A discussion of the limitations 
and possible improvements in the assumptions and the whole approach will be 
presented later. . . • 

Figures 3a, b, and c show the resulting dependence on the. radial distance 
to the center of the nucleus of the coma gas temperature, velocity, and Mach 
number, respectively, for the five coupled gas -dynamic/MCPTM calculations. 

The results have been terminated at the traditional collision zone radius, 
that radial distance to the center of the nucleus which is equal to the local 
mean free path for one collision. Outside this distance, collisions between 
all atoms and molecules are too infrequent to drive the typical adiabatic 
expansion/cooling predicted by pure gas -dynamic models which gradually turns 
internal random kinetic energy (temperature) into larger bulk outflow speeds. 
As discussed in Paper 1, there is a large and gradual transition zone from 
collision-dominated to collisionless molecular flow extending from somewhat 
inside to somewhat outside the traditional collision zone radius. 


The curves presented in Figure 3 contain some significant results in 
terms of both the production of the spatial morphology of the Lyman-a coma and 
also some general properties of cometary comae. In the fairly mild case of a 
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comet with a production rate of 10 s at 1 AU from the sun, the gas 
temperature (after initial cooling) never reaches values greater than about 
30K (Paper 1) . On the other hand, when the heliocentric distance is less than 
0.5 AU, the gas temperatures at some locations within the inner coma can reach 
values ranging from 300 K at 3000km from the nucleus at 0.43 AU to as large as 
900 K at 1000 km from the nucleus at 0.14 AU. Owing to a combination of the 
adiabatic expansion/cooling and the collisional decoupling of the hot H atoms, 
however, the temperatures at the edge of the collision zone radius are always 
smaller than 100 K. 

A heating rate efficiency may be defined as the true local heating rate 
(due to the increasingly insufficient number of collisions by the hot H atoms) 
divided by the rate assuming 100% local thermalization. Figure 5 in Paper 1 
shows that the heating rate efficiency begins to fall from unity at distances 
of about 1/10 of the collision zone radius . Therefore, before the bulk of the 
heavier coma gas can become collisionally decoupled from itself (at the 
collision zone radius) , the photochemical heating rate from the light hot H 
atoms falls and the adiabatic expansion and cooling dominates the flow. So, 
the temperatures at the collision zone radius will generally always be small. 
In terms of their contribution to the velocity distributions of H atoms 
leaving the inner coma, however, the higher temperatures inside the tradi- 
tional collision zone radius are important when considering the large number 
of collisions required to thermalize the atoms compounded by the true direc- 
tional dependence of the collisional mean free path which highly favors radial 
escape (see for example Figure 3 in Paper 1). 



The adopted production rate variation which is asymmetric for equal 
heliocentric distances from pre- to post- perihelion times does have some 
effect of the flow dynamics, although it is only really noticeable for the . 
pair of model results for r - 0.43 AU where the pre-to-post production rate 
ratio is 1.6. At 0.25 AU the ratio is only 1.1 and little difference is 
seen. During the preperihelion times the gas production rate is higher, the 
collision zone radius is larger and the resulting heating rate efficiency 
remains higher for somewhat larger distances from the nucleus. Interestingly 
the same combination discussed above of falling heating efficiency and 
adiabatic expansion and cooling yields a gas temperature that is the same both 
before and after perihelion at the collsion zone radius (which is larger 
before perihelion) . What changes is the bulk outflow speed at the collision 
zone radius which is 1.70 km/s before perihelion and 1.55 km/s afterwards. 

This same process may be at least partially responsible for the different 
outflow speeds inferred from the infrared observations of 1^0 in comet 
P/Halley at comparable heliocentric distances before and after perihelion 
(Larson et al. 1986). 

A more general result present in these model curves is the heliocentric 
distance dependence of the outflow speed of the parent gas molecules . Because 
of the self- decoupling of the coma molecules from themselves near the 
collision zone radius, as discussed above, the velocity distribution of 
molecules outside this radius can be approximately described as the vector sum 
of the bulk radial speed at that radius with a small (<0.3 km/s) initially 
isotropic thermal component. This isotropic component will become less 
isotropic and eventually just a small radial distribution at distances much 
larger than the collision zone radius due to simple geometry. Comparing the 
asymptotic values of the outflow speed at different heliocentric distances we 


find that they do vary approximately as the inverse square root of the 
heliocentric distance as in the Delsemme (1982) law. However, the value at 1 
AU is not 0.58 k m/s but more like 1 to 1.1 km/s. The example gas-dynamic/ 
MCPTM results in Paper 1 at 1 AU are consistent with this value. Recent 
estimates of cometary outflow speeds from various measurements of comet 
P/Halley are also consistent with these larger levels (Larson et al 1986, 
Balsiger et al. 1986). 



IV. The Spatial Morphology of the Lyman-a Coma 

The steady-state inner coma model results, as described in the previous 
section of the paper, give the dependence of the outflow speed and temperature 
of the principal coma gases on the cometocentric distance for five different 
times (or locations) of comet Kohoutek. These five snapshots were used as 
input conditions to characterize the time development of the coma dynamics 
which effects both the source kinematics and the collisional targets at a 
kinetic level, in the full three-dimensional time -dependent MCPTM description 
of the spatial morphology of the hydrogen coma. Time -dependent temperatures 
and outflow speeds at all locations within the coma were generated by inter- 
polating between the five steady-state inner coma MCPTM results. Although the 
mathematical description of the model is given in the accompanying Paper la 
brief qualitative summary of the extended coma MCPTM is given below. 

The full hydrogen coma MCPTM actually contains two major parts. One 
treats the fully time -dependent inner coma in a nucleus -centered coordinate 
system in which displacements from the nucleus are small ( < 10** km) as 
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compared with typical heliocentric distances ( 10 km) . This part contains 
most of the complicated physical processes, such as the collisions, disso- 
ciation times, and excess photodissociation energy, as well as the weighting 
scheme employed in the Monte Carlo simulation. This inner coma portion of the 
full hydrogen coma model is completely separate from the steady- state inner 
coma model discussed in the previous section, although many of the procedures 
employed in the two are similar. The second part explicitly calculates the 
heliocentric trajectories of the H atoms in three dimensions given their 
initial conditions as determined from the inner coma part, and performs the 
sky-plane projection of the atom taking into account the relative geometry of 
the earth, sun and comet at the time of the observation. 



Although this full version of the MCPTM is needed to describe the time 
dependence and three-dimensional trajectories required to explain the 
morphology of the spatially extended hydrogen coma, it should be noted that 
the model actually provides a correct and complete description of the coma 
even for the smaller spatially confined inner regions of the coma. This is an 
important aspect which distinguishes this treatment from other models which 
can only be relevant to either the inner or the extended regions but not 
both. The reason for including this generality is that the inner coma 
conditions are in fact critically important in shaping the spatially extended 
coma. Therefore, it is useful at this point to qualitatively go through a 
step-by-step description of a single pair of hydrogen atom trajectories as 
treated in the full model. 

Many primary parent water molecules are considered to be vaporized from 
the nucleus over a long backup time interval prior to the observation time 
snapshot. The emission time distribution is not skewed; the emission times 
for individual molecules are picked randomly on the backup time interval 
(Combi and Delsemme 1980a) . The variable production rate is included as a 
multiplicative weighting factor for the contribution of both daughter H 
atoms . A forced dissociation time for the H 2 O molecule is calculated 
according to the general time -dependent description for the lifetime, and the 
molecule is displaced from the surface of the nucleus by the appropriate 
distance according to time -dependent coma dynamics. 

At this point the H 2 O molecule is dissociated. .Separate trajectories are 
followed for the OH radical and the primary H atom, as they are ejected from 
the moving center-of-mass of the H 2 O molecule, which is given its appropriate 
local bulk radial speed and a randomized isotropic thermal component from the 
local gas temperature. The ejection speeds of the fragments (H and OH) are 


picked according to the excess energy speed distribution of the photodissocia- 
tion (Table 1). For the OH radical a forced dissociation time is also cal- 
culated using the same method as for the parent ^0, except that the explicit 
heliocentric velocity dependence of the OH lifetime (Schleicher 1983, van 
Dishoeck and Dalgamo 1984) is also included. The OH trajectory is then 
calculated, collisions with the outflowing coma gas molecules are allowed to 
occur, and then finally the radical is dissociated at the appropriate time. 
Another H atom is then emitted from the moving center-of-mass of the OH 
according to the excess energy distribution for that dissociation (Table 2) . 

The trajectories of each of the H atoms (the primary from the H 2 O 
dissociation and the secondary from the OH dissociation) are followed in a 
similar manner to the OH radical within the region where collisions occur. 

A distance of 10** km was adopted as an upper bound distance to separate the 
inner and extended coma portions of the MCPTM. Nearly all of the collisions 
occur at much smaller distances. The H atoms are followed in the model 
initially in the inner. coma portion of the model to calculate the multiple 
collisions which occur. The typical scenario is that at some distance from 
the nucleus which is much smaller than 10^ km the collision path (Eq. 32 in 
Paper 1) becomes infinite; at this point the initial conditions of the atom 
(location, speed and time) are set for the extended coma trajectory 
calculation. This is done by applying a fourth- order Runge-Kutta method to 
the modified classical two-body problem (Eqs. 34-36 in Paper 1). 

For each time step (either an inter-collision time interval or a time 
step in the Runge-Kutta solution of the extended trajectory) the time- 
dependent decay weights for the H atoms are calculated according to Eq. 12 in 
Paper 1. At the observation time the relative projection geometry of the 
earth, sun, comet and atom is calculated and a sky-plane location for the atom 



relative to the comet nucleus is found. For a remote observation, for example 
from Pioneer Venus (Combi, Stewart and Smyth 1986), the location of the earth 
would of course be replaced with the location of Venus. The Lyman-a emission 
rate is calculated for the actual heliocentric velocity and distance of the 
atom according to the adopted solar Lyman -a profile and integrated flux. The 
contribution of the atom is then added to sky-plane grids for density and 
Lyman-a emission having already been defined. Since the trajectories and all 
of the physical processes, at least for the results presented here, are 
symmetric above and below the comet orbit plane each atom is also reflected 
through this plane (z -* -z) and an identical contribution is added at that 
location. Given an arbitrary observational geometry, as is generally the 
case, this provides a direct gain in efficiency of a factor of two. An above 
the orbit plane view of the coma for a fictitous observer located at a large 
distance above and normal to the plane of the comet's orbit, is also saved .in 
the model . This serves as an interesting reference point usually showing the 
most elongated tail. 

The full hydrogen coma MCPTM, as described qualitatively above, and in 
detail in the previous two sections of this paper and in Paper 1, has been 
applied to analyze the two Lyman-a coma images of comet Kohoutek published in 
a paper by Meier et al. (1976). The first of the images was recorded by an 
electrographic camera on board Skylab on December 25, 1973 when the comet was 
0.18 AU from the sun. The second image was recorded by a rocket -borne 
electrographic camera on January 8, 1974 when the comet was 0.43 AU from the 
sun. As discussed in section I, other images recorded from Skylab did not 
produce high quality two-dimensional images and only the upsun and downsun 
profiles were analyzed by Meier et al. with the model of Keller and Meier 


(1976). 



The full isophote contour shapes and their brightness gradients for both 
images were successfully reproduced with the same self-consistent model based 
on a production rate variation determined from the visual light curve (which 
has been reproduced by Delsemme's (1976) standard vaporization model). All 
other physical processes were modeled according to their naturally expected 
heliocentric distance and velocity dependences. Figure 4 shows the comparison 
of the modeled and observed isophote contours for comet Kohoutek on January 8, 
1974. The total hydrogen column densities are small enough (at least up to 
the 50 kR level) so that the coma is clearly optically thin. The small scale 
irregularies in the model contours are generally representative of the 
statistical uncertainties which are unavoidable in this type of analysis. 
However, this poses no problem here since they are much smaller in scale than 
the observed irregularities. 

Figures 5a and b show two versions of the comparison of the modeled 
isophote contours for the observation of December 25. Since the column 
density of hydrogen is large, owing to the large production rate, it is likely 
that the central portion of the coma is somewhat optically thick. It is clear 
from the physical processes occuring in the inner coma and treated by our 
model, that neither a simple spherical radial outflow nor a thermally 
isotropic description (as are typically assumed) will provide the correct H- 
atom velocity distribution which is needed to calculate properly the radiative 
transfer of the Lyman-a photons in the inner coma. The calculation of the 
spherical radiative transfer problem using a real is i tic phase space (space 
density and velocity) distribution for H atoms in the coma (as is included in 
our model) would be an important but difficult problem. 

In order to attempt a first order estimate of the effect of radiative 
transfer on the Lyman-a coma, we have compared in Figure 5 the observations 


both with the optically thin results and with a plane parallel correction to 
the radiative transfer problem (McElroy and Yung, 1975). The optically thin 
result should overestimate the brightness in the inner optically thick region, 
whereas the plane parallel correction should underestimate the brightness. 

This particular radiative transfer correction is for a plane parallel slab of 
hydrogen at a temperature of 5000 K, which corresponds to a most probable 
Maxwell -Boltzmann speed of 9.2 km/s which is actually quite representative of 
the mean speeds in the inner coma. A true spherical radiative transfer calcu- 
lation should fall somewhere intermediate between the plane -parallel and the 
optically thin result. The comparison of the observed image with the two 
modeled images in fact shows that a correction generally intermediate between 
these two extremes would work reasonably well. A realistic radiative transfer 
calculation should account for the true phase space distribution which is 
neither thermally isotropic (like the plane parallel case) nor perfectly 
radial . 

Another complication included in our model results for the December 25, 
1973 image is that for absorption of the cometary Lyman-a emission by the 
hydrogen atoms in the earth's geocorona. On December 25, 1973 the comet's 
geocentric velocity was only 2.29 km/s. Ue have assumed the same transmission 
characteristics as a function of wavelength as used by Meier et al. (1976) in 
their model result for the same image. One improvement enabled by the use of 
the MCPTM method, which provides the actual three -vector velocity of each atom 
trajectory, is that the explicit transmission function was used to attenuate 
the emission contribution for each atom. Meier et al. used an approximate 
form for the assumed width of the velocity distribution of cometary H atoms 
along each column line of sight, given the average doppler shift of the whole 



column from the syndyname model, and then integrated the transmission curve 
convolved with the approximate cometary emission line profile. 

As has already been mentioned the same normalization factor (Qq in Figure 
2) , was used in scaling the variation of the gas production rate to its actual 
time -dependent variation for both image. This normalization factor sets the 

value of the hydrogen atom production rate at 1 AU preperihelion to be 2.7 x 
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10 per second. This implies a water production rate of 1.47 x 10 per 

second at 1 AU pre-perihelion, given the adopted water photodissociation 
branching ratios. This yields hydrogen production rates at the two observa- 
tion times of 9.6 x 10^ at 0.43 AU and 8.4 x 10^® per second at 0.18 AU, as 
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compared with 6.2 x 10 and 5.5 x 10 , respectively, found by Meier et al. 

(1976) using the syndyname model. There is really not much significance in 
the difference in production rate, since most of the difference is due to the 
fact that we adopted a smaller value for the solar Lyman-a flux based on more 
recent solar data. There are, however, a number of significant implications 
which may be drawn .by the comparision of the two very different modeling 
analyses (presented here and by Meier et al,) both of which reproduce the 
shapes of the isophote contours reasonably well. 

First the fundamental differences between the two models should again be 
stated. The syndyname models used by Meier et al. consider the source of 
hydrogen to be a point centered on the nucleus. For the purpose of analyzing 
these wide field images, this in itself is a reasonable simplification since 
the scale lengths for the source dissociations of l^O and OH are not 
resolvable on the scale of the data. Their outflow speed distribution on the 
other hand is parameterized and constructed by the weighted sum of three 
Maxwell -Boltzmann distributions which are, for a single image, assumed to be 
constant in time. For the MCPTM results presented here the effective outflow 
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speed distribution for H atoms leaving the inner coma is not a fitting para- 
meter but represents the natural self-consistent outcome of photodissociations 
and multiple collisions with a water dominated coma. Furthermore, this 
effective speed distribution in the MCPTM is inherently a time dependent 
quantity (unlike the syndyname model) since the collision rates depend in a 
complex way on the production rate, the photochemical lifetimes and the coma 
temperatures, all of which are explicitly time dependent. 

A rough quantitative comparison can be made for the speed, distribution 
functions of hydrogen atoms leaving the inner coma of the comet. It is, after 
all, this speed distribution, which when interacting with the relative forces 
of solar gravity and radiation pressure between different atoms, that creates 
the observed shapes of the isophote contours. Figure 6 compares the sums of 
the three Maxwell -Boltzmann distributions as determined by the fitting of the 
syndyname models by Meier et.al. to the observed images, with the effective 
time-averaged outflow speed distribution function as calculated in the 
extended coma MCPTM results. This distribution was calculated in the MCPTM by 
determining the outflow speed of each hydrogen atom at the time the atom left 
the inner coma portion of the model (always < 10^ km). The outflow speed 
distribution was binned in 1 km/s intervals and each atom's density 
contribution to the coma as a whole was also accumulated in the correct 
velocity bin. This was found to be the most directly comparable quantity to 
compare the MCPTM results with the parameterized distributions of Meier et al. 

When comparing the syndyname results for the two different images, Meier 
et al. noted that a signf leant population of low speed ( 4 km/s component) H 
atoms were needed in order to fit the image at 0.18 AU, but that only a small 
(and possibly negligible) amount was needed in order to explain the image at 
0.43 AU. As the MCPTM results demonstrated, they correctly assessed the fact 



that more H atoms are produced deep in the dense part of the coma where they 
would be thermal ized. High speed components ( > 20 km/s) are also evident in 
both sets of modeled distributions, although since we have adopted a smaller 
solar Lyman-a flux we have obtained a reasonable fit with essentially no atoms 
above 32 km/s. Both sets of model results also peak at around 8 km/s. The 
MCPTM distributions are however far more irregular than the smooth Mawxell- 
Boltzmann distributions. Although, it should be noted that Meier et al. 
explicitly mentioned there was no implication of a unique fit for their 
modeled distributions to the image data. 

As discussed in section II and Paper 1, the assumption of elastic hard- 
sphere collisions between atoms and molecules with cross sectional diameters 
of order of a few Angtrom units is rather simplistic. Actual interactions are 
rather complex, having long range but weak attractive forces and short range 
repulsive forces (Johnson, 1982). In addition, the collision between an atom 
and a molecule is likely to be somewhat inelastic (at least in the general 
case) owing to the possibility of exciting internal rotational and/or 
vibrational states of the molecule. In the absence of measured scattering 
cross sections between H atoms and H 2 O molecules in the appropriate energy 
range (0.1 to 2 eV) , some general statements can be made regarding the 
assumption of elastic hard-sphere collisions and its effect on our model. 

For the purposes of our calculation, if inelastic processes were very 
important then more energy could be transfered per collision. Any excess 
energy transfer from the H atom to a water molecule, above the elastic portion 
which goes initially into tranlational energy, would be absorbed by the water 
molecule in the form of a rotational (or vibrational) excitation. In most of 
the coma, which is optically thin in the infrared and where the intercollsion 
times between water molecules is very long compared with radiative lifetimes, 
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this excess tranferred energy would simply radiate away. The major effect 
would be on the more rapid thermalization of H atoms. This would result in a 
larger number of very low speed atoms than are evident (see Figure 6) . 

Inelastic collision cross sections are generally smaller than elastic 
cross sections (Bernstein, 1979). Such a generalization can also be inferred 
by the fairly long relaxation times required for polyatomic molecules to reach 
complete thermal equilibrium in their rotational and vibrational state 
distributions. For example up to 50 collisions are required for H 2 to reach 
rotational equilibrium, whereas thousands of collisions are required for 
triatomic molecules to reach vibrational equilibrium (Hirschfelder , Curtiss, 
and Bird, 1954). In addition, elastic cross sections are determined by using 
simple elastic models to interpret measured macroscopic quantities such as 
diffusion and viscosity, and thus likely incorporate some kind of average 
inelastic component at least in some implicit way.- Therefore, our use of 
hard- sphere collisions to infer such macroscopic quantities as the energy 
transfer and velocity distribution of H atoms should at least be self- 
consistent. Similar reasoning is also employed in hydrodynamic models 
(Marconi and Mendis, 1983, 1984, and Gombosi et al. 1986) which use momentum 
transfer rates based on the same type of hard-sphere elastic collisions.. 

Another important implication of the reproducibility of the isophote 
contours by the extended coma MCTPM is that we have used the excess speed 
distribution for H atoms dissociated for OH radicals according to the ab 
initio calculations of van Dishoeck and Dalgamo (1984) who predict high lying 
predissociation states, with large excess energies yielding fast H atoms, 
having speeds of 17 to 26 km/s. If we ignore these newly identified states 
and consider all of the H atoms produced from OH only at 8 km/s (as is a 
common practice) the shapes of the modeled isophote contours deviate from the 
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observed shapes significantly. This result and the previously mentioned 
observation of the emission at 6300 A from 0(^D) atoms at nearly 10^ km from 
comet P/Halley by Roesler et al. (1987) would seem to verify the importance of 
these new predissociation transitions. 

The fact that our model can directly account for the speed distribution 
of H atoms leaving the inner coma which successfully reproduces the Lyman-a 
isophotes has some important implications regarding the proper modeling of 
hydrogen in other situations. Keller and Meier (1980) also used the syndyname 
model in order to reproduce the observed Lyman-a images of comet West (1976 
VI) , but they required the presence of an outburst in the production of H at 
some past time. The comet West observation was, like Kohoutek, of a very 
productive comet (actually more productive) at a small heliocentric distance 
(0.4 AU) where collisions should be expected to be important in determining 
the exit speeds of H atoms from the inner coma. However, for comets either 
less productive than Kohoutek at small heliocentric distances, or as 
productive but at larger heliocentric distances, the role of collisions in 
shaping the effective exit speed distribution of H atoms from the inner coma 
would be minimal. 

The effective speed distributions for H atoms in these essentially 
"collisionless" cases would be set by the proper vector sums of several 
components: (1) the speed from the excess energy of all photodissociations, 

(2) the outflow speed of the coma, and (3) the temperature of the coma gas. 

An example of such a "collisionless” speed distribution can be illustrated by 
the dashed histograms in Figure 6 which show the initial speed distributions 
of the H atoms in the Kohoutek models upon their initial production but before 
any collisions. Example models have been presented earlier (Combi and Smyth 
1987) for the case of the rocket image of Kohoutek at 0.43 AU, where the 


production rate was reduced by factors of 3 and 10, thereby reducing signifi- 
cantly the number and importance of collisions. The spatial character of the 
cloud also changed markedly for these cases. Previous MCPTM results published 
for the Pioneer Venus Orbiter UVS observations of comet P/Giacobini-Zinner 
made during the ICE encounter were made using an (appropriately) collisionless 

version of the model which should certainly be valid for this comet at 1.03 AU 
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from the sun where the water production rate was only 2.3 x 10 molecules per 
second (Combi, Stewart and Smyth 1986). 
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V. Discussion 

The application of the MCPTM, described in a general mathematical form in 
the accompanying paper (Paper 1) , has been applied in this paper to explain 
the observed spatial morphology of the Lyman-a coma. The two electrographic 
camera images of comet Kohoutek (1973 XII) were chosen for this particular 
study. Although these images had been reasonably well reproduced in the 
original observation paper by Meier et al. (1976), the physical basis of their 
sums of two or three Maxwell -Boltzmann distributions for the effective H atom 
radial outflow speeds remained unquantified. In this paper the MCPTM accounts 
explicitly for the time -dependent physical state of the inner coma source, the 
proper isotropic (vectorial) ejection of daughter OH radicals and H atoms, the 
collisional thermal izat ion of the H atoms (which is generally only partial) , 
and the explicit calculation of the trajectories of the appropriately produced 
H atoms under the influence of solar radiation pressure and gravity in three 
dimensions . 

While a comparison of the multi -component Maxwellian distributions of 
Meier et al. with the effective outflow speed distributions of H atoms from 
the inner coma as modeled by the MCPTM shows clear differences in the detailed 
appearance, there is enough broad similarity to see why both models can 
produce a reasonable fit to the data. The success of the physically-based 
MCPTM, however, implies that the models of Keller and Meier (1976) cannot be 
applied indiscriminately to the hydrogen observations of comets other than 
very productive comets at small heliocentric distances. In fact, the reason 
the multi -component Maxwellian models are successful is simply due to the 
observational selection effect that only very productive comets at small 
heliocentric distances have been observed in a manner which yield the wide 
field images available for comets Kohoutek and West. Furthermore, it is only 



for comets at small heliocentric distances that the radiation pressure is 
large enough so that appropriate modeling is particularly sensitive to the 
effective outflow speed distribution of H atoms from the inner coma. 

It is important to stress that we have been able to produce an effective 
time -dependent outflow speed distribution of H atoms leaving the inner coma 
which successfully reproduces the observed Lyman-a isophote gradient and 
shapes, using an inherently physically-based self-consistent model. Although 
many processes were included in the model, and best estimates had to be made 
for a number of physical parameters, there are few truely free fitting 
parameters in the entire procedure. Sensitivity studies of the important 
parameters would certainly represent interesting, although somewhat compu- 
tationally intensive, future undertakings. Apart from this there are two 
important areas for future improvements and generalizations which can be made 
to the model and the modeling process. 

In the details of the collision calculation, an obvious improvement would 
be to include the exact expression for the collision path integral (Eq. 31, in 
Paper 1). Although for the first several collisions, where an atom's speed is 
large compared with the outflow speed and where the bulk of the heating takes 
place, the approximate expression (Eq. 32 in Paper 1) is clearly a reasonable 
choice. After a number of collisions, the atoms are more or less thermalized, 
and the exact expression for the collision path is more appropriate. However, 
it will not make much difference to the final atom kinematics if the atoms are 
collided a few too many times once they are already thermalized. Another 
improvement in the area of the collision process is to carefully consider the 
applicability of the oversimplified hard-sphere scattering which has always 
been used (and is here) in order to understand the role of collisions between 
atoms, radicals and molecules. 



The more fundamental direct improvements are to be made to the physical 
inner coma description. The possible effect of dust, plasma, and the minor 
heating mechanisms should be addressed more carefully. The collisional 
decoupling of the hot H atoms from the heavier coma gas molecules represents a 
first order correction to the general problem of the breakdown of fluid 
mechanics in describing not only the correct photochemical heating but more 
generally the physics of the outflow of the cometary coma. A better treatment 
of the problem would be accomplished by including separate MCPTM calculations 
for each of the important species, namely H, OH, 0, ^0, CO and CO 2 . In this 
way not only can a complete description for observable species be made, but 
also, a basic study at the kinetic theory level of the breakdown of the 
applicability of fluid mechanics could be studied. 
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TABLE 3 


Coupled Gas - Dynamic/MCPTM 


Heliocentric 

Distance 

(AU) 


Water Production 


Rate 

-.(s' 1 ) 


Initial Gas e 
Temperature 

no 


0.434 a 
0. 25 a 
0 . 142 b 
0.25° 
0.434 c 


8.8 x 10 29 
2.7 x 10 30 

6.4 x 10 30 

2.5 x 10 30 

5.5 x 10 29 


208 

216 

225 

216 

208 


a pr e - per ihe 1 ion 

^perihelion 

c post-perihelion 

d Based on Q fl - 1.47 x 10 29 s* 1 at 
e Extrapolated from Gombosi et al. 


1 AU pre-perihelion in Figure 2. 
(1985) 


Figure Captions 


Figure 1. Solar Lyman-a Line Profile. The full solar disc Lyman-a line 
profile as adopted in the MCPTM was constructed by Lemaire et al. (1978) from 
high resolution disc -centered and limb spectra. The geocoronal absorption at 
line center was filled- in to approximate the flux as seen by a hydrogen atom 
in the cometary coma. 

Figure 2. Variation of the Relative Gas Production Rate in Comet 
Kohoutek. The heliocentric distance dependence of the relative gas production 
rate as adopted from the visual light curve of comet Kohoutek which was com- 
piled by Delsemme (1976). The variation is normalized to the production rate 
at 1 AU pre -perihelion. The circles show the pre -perihelion values, while the 
stars show the post-perihelion values. 

Figure 3. Five Coupled Gas -dynamic/MCPTM Coma Results for Comet 
Kohoutek. Shown are the temperature, outflow speed, and Mach number of the 
coma gas as a function of radial distance to the nucleus. The dashed curves 
show the model results for the pre-perihelion heliocentric distances indi- 
cated, whereas the solid lines give the post-perihelion results. The coma 
model parameters are given in Table 3 and all model curves are terminated at 
the boundary of the collision zone radius. 

Figure 4. MCPTM/Data Comparison of the Lyman-a Coma of Comet Kohoutek at 
0.43 AU. The isophote contours for the data (solid lines) as observed by 
Meier et al. (1976) are compared with the fully time -dependent 3D self- 
consistent MCPTM result (dashed lines) . The contours are in units of 
kiloRayleighs and the tick marks on the circumscribed box are separated by 



Figure 5. MCPTM/Data Comparisons of the Lyman-a Coma of Comet Kohoutek 
at 0.18 AU. The isophote contours for the data (solid lines) as observed by 
Meier et al. (1976) are compared with the fully time -dependent 3D self- 
consistent MCPTM result (dashed lines). Part (a) shows the optically thin 
result, and part (b) shows the result of a simple plane parallel radiative 
transfer correction (T-5000K) to the thin case. A true spherical radiative 
transfer calculation which takes the real phase space distribution (density 
and velocity) into account should fall intermediately between these two 
models. The contours are in units of kiloRayleighs and the tick marks on the 
circumscribed box are separated by 4.0 x 10^ km. 

Figure 6. A Comparison of the Velocity Distribution Functions for H 
atoms Leaving the Inner Coma in the MCPTM and Syndyname Model . The smooth 
solid curves show the sums of the three -component Maxwellians of the syndyname 
models of Meier et al. (1976). The solid histograms give the MCPTM distri- 
bution functions for H atoms leaving the inner collisional coma and transition 
zone. The dashed histograms give the MCPTM distribution functions for the H 
atoms upon initial production at dissociation of the parent H 2 O molecules and 
OH radicals but before collisions. Part (a) gives the results for comet 
Kohoutek at 0.18 AU, and part (b) gives the results at 0.43 AU. 
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